diff --git a/Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc b/Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc index 5bfa958662f29..22e4050ef707a 100644 --- a/Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc +++ b/Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc @@ -55,8 +55,9 @@ #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h" +#include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h" -#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h" +#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h" #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h" #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" #include "TrackingTools/TransientTrack/interface/TransientTrack.h" @@ -120,11 +121,6 @@ PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfi theTrackClusterizer_ = std::make_unique(iConfig.getParameter("TkClusParameters") .getParameter("TkGapClusParameters")); - } else if (clusteringAlgorithm == "DA") { - theTrackClusterizer_ = - std::make_unique(iConfig.getParameter("TkClusParameters") - .getParameter("TkDAClusParameters")); - // provide the vectorized version of the clusterizer, if supported by the build } else if (clusteringAlgorithm == "DA_vect") { theTrackClusterizer_ = std::make_unique(iConfig.getParameter("TkClusParameters") @@ -3700,7 +3696,7 @@ void PrimaryVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& d // track filtering edm::ParameterSetDescription psd0; TrackFilterForPVFinding::fillPSetDescription(psd0); - psd0.add("numTracksThreshold", 0); // HI only + HITrackFilterForPVFinding::fillPSetDescription(psd0); // HI only desc.add("TkFilterParameters", psd0); // PV Clusterization @@ -3708,7 +3704,7 @@ void PrimaryVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& d edm::ParameterSetDescription psd0; { edm::ParameterSetDescription psd1; - DAClusterizerInZ_vect::fillPSetDescription(psd1); + DAClusterizerInZT_vect::fillPSetDescription(psd1); psd0.add("TkDAClusParameters", psd1); edm::ParameterSetDescription psd2; diff --git a/Alignment/OfflineValidation/python/TkAlAllInOneTool/PV_cfg.py b/Alignment/OfflineValidation/python/TkAlAllInOneTool/PV_cfg.py index 071ba09157af5..33ca243a73a24 100644 --- a/Alignment/OfflineValidation/python/TkAlAllInOneTool/PV_cfg.py +++ b/Alignment/OfflineValidation/python/TkAlAllInOneTool/PV_cfg.py @@ -187,7 +187,7 @@ ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Alignment/OfflineValidation/test/PVValidation_TEMPL_cfg.py b/Alignment/OfflineValidation/test/PVValidation_TEMPL_cfg.py index 1b5482137751a..968ad2c66aac6 100644 --- a/Alignment/OfflineValidation/test/PVValidation_TEMPL_cfg.py +++ b/Alignment/OfflineValidation/test/PVValidation_TEMPL_cfg.py @@ -216,7 +216,7 @@ def customiseKinksAndBows(process): ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Alignment/OfflineValidation/test/PVValidation_T_cfg.py b/Alignment/OfflineValidation/test/PVValidation_T_cfg.py index 866e9413c8235..3349d4acc1879 100644 --- a/Alignment/OfflineValidation/test/PVValidation_T_cfg.py +++ b/Alignment/OfflineValidation/test/PVValidation_T_cfg.py @@ -279,7 +279,7 @@ def customiseKinksAndBows(process): ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Alignment/OfflineValidation/test/testG4Refitter_cfg.py b/Alignment/OfflineValidation/test/testG4Refitter_cfg.py index 050c40634534e..ab79fb97960a3 100644 --- a/Alignment/OfflineValidation/test/testG4Refitter_cfg.py +++ b/Alignment/OfflineValidation/test/testG4Refitter_cfg.py @@ -130,7 +130,7 @@ ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() #################################################################### diff --git a/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py b/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py index 8ec25c822b0ca..d4ef8bc13880c 100644 --- a/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py +++ b/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py @@ -285,7 +285,7 @@ class RefitType(Enum): ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Alignment/OfflineValidation/test/test_all_Phase2_cfg.py b/Alignment/OfflineValidation/test/test_all_Phase2_cfg.py index cca0e44db5153..1c9a437574dc1 100644 --- a/Alignment/OfflineValidation/test/test_all_Phase2_cfg.py +++ b/Alignment/OfflineValidation/test/test_all_Phase2_cfg.py @@ -258,7 +258,7 @@ class RefitType(Enum): ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Alignment/OfflineValidation/test/test_all_cfg.py b/Alignment/OfflineValidation/test/test_all_cfg.py index 024b97bd3d9d7..c6bc6945fdbe5 100644 --- a/Alignment/OfflineValidation/test/test_all_cfg.py +++ b/Alignment/OfflineValidation/test/test_all_cfg.py @@ -256,7 +256,7 @@ class RefitType(Enum): ) ## MM 04.05.2017 (use settings as in: https://github.com/cms-sw/cmssw/pull/18330) -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import DA_vectParameters DAClusterizationParams = DA_vectParameters.clone() GapClusterizationParams = cms.PSet(algorithm = cms.string('gap'), diff --git a/Configuration/ProcessModifiers/python/vertex4DTrackSelMVA_cff.py b/Configuration/ProcessModifiers/python/vertex4DTrackSelMVA_cff.py new file mode 100644 index 0000000000000..4983c6b1fe0c3 --- /dev/null +++ b/Configuration/ProcessModifiers/python/vertex4DTrackSelMVA_cff.py @@ -0,0 +1,6 @@ +import FWCore.ParameterSet.Config as cms + +# Modifier to enable the use of the MVA selection on +# tracks for the 4D vertex reco + +vertex4DTrackSelMVA = cms.Modifier() diff --git a/RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc b/RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc index 1ae9f060d8fa6..d1144ca7c2ec1 100644 --- a/RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc +++ b/RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc @@ -49,11 +49,18 @@ class TOFPIDProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT> tofkToken_; edm::EDGetTokenT> tofpToken_; edm::EDGetTokenT vtxsToken_; - double vtxMaxSigmaT_; - double maxDz_; - double maxDtSignificance_; - double minProbHeavy_; - double fixedT0Error_; + edm::EDGetTokenT> trackMTDTimeQualityToken_; + const double vtxMaxSigmaT_; + const double maxDz_; + const double maxDtSignificance_; + const double minProbHeavy_; + const double fixedT0Error_; + const double probPion_; + const double probKaon_; + const double probProton_; + const double minTrackTimeQuality_; + const bool MVASel_; + const bool vertexReassignment_; }; TOFPIDProducer::TOFPIDProducer(const ParameterSet& iConfig) @@ -65,11 +72,19 @@ TOFPIDProducer::TOFPIDProducer(const ParameterSet& iConfig) tofkToken_(consumes>(iConfig.getParameter("tofkSrc"))), tofpToken_(consumes>(iConfig.getParameter("tofpSrc"))), vtxsToken_(consumes(iConfig.getParameter("vtxsSrc"))), + trackMTDTimeQualityToken_( + consumes>(iConfig.getParameter("trackMTDTimeQualityVMapTag"))), vtxMaxSigmaT_(iConfig.getParameter("vtxMaxSigmaT")), maxDz_(iConfig.getParameter("maxDz")), maxDtSignificance_(iConfig.getParameter("maxDtSignificance")), minProbHeavy_(iConfig.getParameter("minProbHeavy")), - fixedT0Error_(iConfig.getParameter("fixedT0Error")) { + fixedT0Error_(iConfig.getParameter("fixedT0Error")), + probPion_(iConfig.getParameter("probPion")), + probKaon_(iConfig.getParameter("probKaon")), + probProton_(iConfig.getParameter("probProton")), + minTrackTimeQuality_(iConfig.getParameter("minTrackTimeQuality")), + MVASel_(iConfig.getParameter("MVASel")), + vertexReassignment_(iConfig.getParameter("vertexReassignment")) { produces>(t0Name); produces>(sigmat0Name); produces>(t0safeName); @@ -97,6 +112,8 @@ void TOFPIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio ->setComment("Input ValueMap for track tof as proton"); desc.add("vtxsSrc", edm::InputTag("unsortedOfflinePrimaryVertices4DwithPID")) ->setComment("Input primary vertex collection"); + desc.add("trackMTDTimeQualityVMapTag", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA")) + ->setComment("Track MVA quality value"); desc.add("vtxMaxSigmaT", 0.025) ->setComment("Maximum primary vertex time uncertainty for use in particle id [ns]"); desc.add("maxDz", 0.1) @@ -107,6 +124,12 @@ void TOFPIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio desc.add("minProbHeavy", 0.75) ->setComment("Minimum probability for a particle to be a kaon or proton before reassigning the timestamp"); desc.add("fixedT0Error", 0.)->setComment("Use a fixed T0 uncertainty [ns]"); + desc.add("probPion", 1.)->setComment("A priori probability pions"); + desc.add("probKaon", 1.)->setComment("A priori probability kaons"); + desc.add("probProton", 1.)->setComment("A priori probability for protons"); + desc.add("minTrackTimeQuality", 0.8)->setComment("Minimum MVA Quality selection on tracks"); + desc.add("MVASel", false)->setComment("Use MVA Quality selection"); + desc.add("vertexReassignment", true)->setComment("Track-vertex reassignment"); descriptions.add("tofPIDProducer", desc); } @@ -142,6 +165,8 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) { const auto& vtxs = ev.get(vtxsToken_); + const auto& trackMVAQualIn = ev.get(trackMTDTimeQualityToken_); + //output value maps (PID probabilities and recalculated time at beamline) std::vector t0OutRaw; std::vector sigmat0OutRaw; @@ -165,7 +190,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) { float prob_k = -1.; float prob_p = -1.; - if (sigmat0 > 0.) { + float trackMVAQual = trackMVAQualIn[trackref]; + + if (sigmat0 > 0. && (!MVASel_ || (MVASel_ && trackMVAQual >= minTrackTimeQuality_))) { double rsigmazsq = 1. / track.dzError() / track.dzError(); double rsigmat = 1. / sigmatmtd; @@ -239,7 +266,12 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) { double chisqmin_k = std::numeric_limits::max(); double chisqmin_p = std::numeric_limits::max(); //loop through vertices and check for better matches - for (const reco::Vertex& vtx : vtxs) { + for (unsigned int ivtx = 0; ivtx < vtxs.size(); ++ivtx) { + const reco::Vertex& vtx = vtxs[ivtx]; + if (!vertexReassignment_) { + if (ivtx != (unsigned int)vtxidx) + continue; + } if (!(vtx.tError() > 0. && vtx.tError() < vtxMaxSigmaT_)) { continue; } @@ -283,9 +315,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) { //compute PID probabilities //*TODO* deal with heavier nucleons and/or BSM case here? - double rawprob_pi = exp(-0.5 * chisqmin_pi); - double rawprob_k = exp(-0.5 * chisqmin_k); - double rawprob_p = exp(-0.5 * chisqmin_p); + double rawprob_pi = probPion_ * exp(-0.5 * chisqmin_pi); + double rawprob_k = probKaon_ * exp(-0.5 * chisqmin_k); + double rawprob_p = probProton_ * exp(-0.5 * chisqmin_p); double normprob = 1. / (rawprob_pi + rawprob_k + rawprob_p); diff --git a/RecoVertex/Configuration/python/RecoVertex_EventContent_cff.py b/RecoVertex/Configuration/python/RecoVertex_EventContent_cff.py index 00327cd5ec7a3..fc0802a044dd1 100644 --- a/RecoVertex/Configuration/python/RecoVertex_EventContent_cff.py +++ b/RecoVertex/Configuration/python/RecoVertex_EventContent_cff.py @@ -17,9 +17,7 @@ 'keep *_offlinePrimaryVertices4DWithBS__*', 'keep *_trackTimeValueMapProducer_*_*' ] -_phase2_tktiming_layer_RecoVertexEventContent = [ 'keep *_offlinePrimaryVertices4DnoPID__*', - 'keep *_offlinePrimaryVertices4DnoPIDWithBS__*', - 'keep *_tofPID_*_*'] +_phase2_tktiming_layer_RecoVertexEventContent = [ 'keep *_tofPID_*_*'] phase2_timing.toModify( RecoVertexAOD, outputCommands = RecoVertexAOD.outputCommands + _phase2_tktiming_RecoVertexEventContent) phase2_timing_layer.toModify( RecoVertexAOD, diff --git a/RecoVertex/Configuration/python/RecoVertex_cff.py b/RecoVertex/Configuration/python/RecoVertex_cff.py index 2bbba3b2aab00..dd2b3f8927818 100644 --- a/RecoVertex/Configuration/python/RecoVertex_cff.py +++ b/RecoVertex/Configuration/python/RecoVertex_cff.py @@ -45,15 +45,11 @@ from RecoVertex.Configuration.RecoVertex_phase2_timing_cff import (tpClusterProducer , quickTrackAssociatorByHits , trackTimeValueMapProducer , - unsortedOfflinePrimaryVertices4DnoPID , - trackWithVertexRefSelectorBeforeSorting4DnoPID , - trackRefsForJetsBeforeSorting4DnoPID , - offlinePrimaryVertices4DnoPID , - offlinePrimaryVertices4DnoPIDWithBS, unsortedOfflinePrimaryVertices4DwithPID , offlinePrimaryVertices4DwithPID , offlinePrimaryVertices4DwithPIDWithBS, tofPID, + tofPID3D, tofPID4DnoPID, unsortedOfflinePrimaryVertices4D, trackWithVertexRefSelectorBeforeSorting4D, @@ -73,11 +69,7 @@ ) _phase2_tktiming_layer_vertexrecoTask = cms.Task( _phase2_tktiming_vertexrecoTask.copy() , - unsortedOfflinePrimaryVertices4DnoPID , - trackWithVertexRefSelectorBeforeSorting4DnoPID , - trackRefsForJetsBeforeSorting4DnoPID , - offlinePrimaryVertices4DnoPID , - offlinePrimaryVertices4DnoPIDWithBS, + tofPID3D, tofPID, tofPID4DnoPID, ) @@ -93,3 +85,6 @@ phase2_timing_layer.toModify(offlinePrimaryVertices4D, vertices = "unsortedOfflinePrimaryVertices4D", particles = "trackRefsForJetsBeforeSorting4D") phase2_timing_layer.toModify(offlinePrimaryVertices4DWithBS, vertices = "unsortedOfflinePrimaryVertices4D:WithBS", particles = "trackRefsForJetsBeforeSorting4D") +from Configuration.ProcessModifiers.vertex4DTrackSelMVA_cff import vertex4DTrackSelMVA +vertex4DTrackSelMVA.toModify(unsortedOfflinePrimaryVertices4D, useMVACut = True) +vertex4DTrackSelMVA.toModify(unsortedOfflinePrimaryVertices4DwithPID, useMVACut = True) diff --git a/RecoVertex/Configuration/python/RecoVertex_phase2_timing_cff.py b/RecoVertex/Configuration/python/RecoVertex_phase2_timing_cff.py index 1ca0af44b269e..37f196eabaa1b 100644 --- a/RecoVertex/Configuration/python/RecoVertex_phase2_timing_cff.py +++ b/RecoVertex/Configuration/python/RecoVertex_phase2_timing_cff.py @@ -1,12 +1,19 @@ import FWCore.ParameterSet.Config as cms from RecoVertex.Configuration.RecoVertex_cff import unsortedOfflinePrimaryVertices, trackWithVertexRefSelector, trackRefsForJets, sortedPrimaryVertices, offlinePrimaryVertices, offlinePrimaryVerticesWithBS,vertexrecoTask -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA2D_vectParameters - unsortedOfflinePrimaryVertices4D = unsortedOfflinePrimaryVertices.clone( - TkClusParameters = DA2D_vectParameters, + TkClusParameters = dict( + algorithm = "DA2D_vect", + TkDAClusParameters = dict( + Tmin = 4.0, + Tpurge = 4.0, + Tstop = 2.0 + ), + ), TrackTimesLabel = cms.InputTag("trackTimeValueMapProducer","generalTracksConfigurableFlatResolutionModel"), TrackTimeResosLabel = cms.InputTag("trackTimeValueMapProducer","generalTracksConfigurableFlatResolutionModelResolution"), + vertexCollections = {0: dict(vertexTimeParameters = cms.PSet( algorithm = cms.string('legacy4D'))), + 1: dict(vertexTimeParameters = cms.PSet( algorithm = cms.string('legacy4D')))} ) trackWithVertexRefSelectorBeforeSorting4D = trackWithVertexRefSelector.clone( vertexTag = "unsortedOfflinePrimaryVertices4D", @@ -27,28 +34,6 @@ vertices = "unsortedOfflinePrimaryVertices4D:WithBS" ) -unsortedOfflinePrimaryVertices4DnoPID = unsortedOfflinePrimaryVertices4D.clone( - TrackTimesLabel = "trackExtenderWithMTD:generalTrackt0", - TrackTimeResosLabel = "trackExtenderWithMTD:generalTracksigmat0" -) -trackWithVertexRefSelectorBeforeSorting4DnoPID = trackWithVertexRefSelector.clone( - vertexTag = "unsortedOfflinePrimaryVertices4DnoPID", - ptMax = 9e99, - ptErrorCut = 9e99 -) -trackRefsForJetsBeforeSorting4DnoPID = trackRefsForJets.clone( - src = "trackWithVertexRefSelectorBeforeSorting4DnoPID" -) -offlinePrimaryVertices4DnoPID = offlinePrimaryVertices4D.clone( - vertices = "unsortedOfflinePrimaryVertices4DnoPID", - particles = "trackRefsForJetsBeforeSorting4DnoPID", - trackTimeTag = "trackExtenderWithMTD:generalTrackt0", - trackTimeResoTag = "trackExtenderWithMTD:generalTracksigmat0" -) -offlinePrimaryVertices4DnoPIDWithBS=offlinePrimaryVertices4DnoPID.clone( - vertices = "unsortedOfflinePrimaryVertices4DnoPID:WithBS" -) - unsortedOfflinePrimaryVertices4DwithPID = unsortedOfflinePrimaryVertices4D.clone( TrackTimesLabel = "tofPID4DnoPID:t0safe", TrackTimeResosLabel = "tofPID4DnoPID:sigmat0safe" @@ -76,9 +61,15 @@ from SimTracker.TrackAssociation.trackTimeValueMapProducer_cfi import trackTimeValueMapProducer from RecoMTD.TimingIDTools.tofPIDProducer_cfi import tofPIDProducer -tofPID4DnoPID=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices4DnoPID') +tofPID4DnoPID=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices') tofPID=tofPIDProducer.clone() +tofPID3D=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices') from Configuration.Eras.Modifier_phase2_timing_layer_cff import phase2_timing_layer -phase2_timing_layer.toModify(tofPID, vtxsSrc='unsortedOfflinePrimaryVertices4D') +phase2_timing_layer.toModify(tofPID, vtxsSrc='unsortedOfflinePrimaryVertices4D', vertexReassignment=False) +phase2_timing_layer.toModify(tofPID3D, vertexReassignment=False) +phase2_timing_layer.toModify(unsortedOfflinePrimaryVertices, + vertexCollections = {0: dict(vertexTimeParameters = cms.PSet( algorithm = cms.string('fromTracksPID'))), + 1: dict(vertexTimeParameters = cms.PSet( algorithm = cms.string('fromTracksPID')))} +) diff --git a/RecoVertex/PrimaryVertexProducer/interface/AdaptiveChisquarePrimaryVertexFitter.h b/RecoVertex/PrimaryVertexProducer/interface/AdaptiveChisquarePrimaryVertexFitter.h new file mode 100644 index 0000000000000..612d6a8db6b89 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/AdaptiveChisquarePrimaryVertexFitter.h @@ -0,0 +1,98 @@ +#ifndef RecoVertex_PrimaryVertexProducer_AdaptiveChisquarePrimaryVertexFitter_h +#define RecoVertex_PrimaryVertexProducer_AdaptiveChisquarePrimaryVertexFitter_h + +/**\class AdaptiveChisquarePrimaryVertexFitter + + Description: simultaneous chisquared fit of primary vertices + +*/ +#include + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h" + +class AdaptiveChisquarePrimaryVertexFitter : public PrimaryVertexFitterBase { +public: + AdaptiveChisquarePrimaryVertexFitter(double chicutoff = 2.5, + double zcutoff = 1.0, + double mintrkweight = 0.4, + bool multivertexfit = false); + ~AdaptiveChisquarePrimaryVertexFitter() override = default; + + std::vector fit(const std::vector &, + const std::vector &, + const reco::BeamSpot &, + const bool) override; + + using Error3 = ROOT::Math::SMatrix; + +protected: + void verify() { // DEBUG only + unsigned int nt = trackinfo_.size(); + unsigned int nv = xv_.size(); + assert((yv_.size() == nv) && "yv size"); + assert((zv_.size() == nv) && "zv size"); + assert((tkfirstv_.size() == (nv + 1)) && "tkfirstv size"); + assert((tkmap_.size() == tkweight_.size()) && "tkmapsize <> tkweightssize"); + for (unsigned int k = 0; k < nv; k++) { + assert((tkfirstv_[k] < tkweight_.size()) && "tkfirst[k]"); + assert((tkfirstv_[k + 1] <= tkweight_.size()) && "tkfirst[k+1]"); + assert((tkfirstv_[k] <= tkfirstv_[k + 1]) && "tkfirst[k/k+1]"); + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; k++) { + assert((j < tkmap_.size()) && "illegal tkfirst entry"); + unsigned int i = tkmap_[j]; + assert((i < nt) && "illegal tkmap entry"); + assert((tkweight_[i] >= 0) && "negative tkweight or nan"); + assert((tkweight_[i] <= 1) && "tkweight > 1 or nan"); + } + } + }; + + struct TrackInfo { + double S11, S22, S12; // inverse of the covariance (sub-)matrix + Error3 C; // H^T S H + double g[3]; + double H1[3], H2[3]; + double b1, b2; + double zpca, dzError; + }; + + std::vector vertices(const std::vector &, + const std::vector &, + const reco::BeamSpot &, + const bool); + TransientVertex refit(const TransientVertex &, const reco::BeamSpot &, const bool); + double track_in_vertex_chsq(const TrackInfo &, const double, const double, const double); + void fill_trackinfo(const std::vector &, const reco::BeamSpot &); + void fill_weights(const reco::BeamSpot &, const double beta = 1.); + TransientVertex get_TransientVertex(const unsigned int, + const std::vector> &, + const std::vector &, + const float, + const reco::BeamSpot &); + Error3 get_inverse_beam_covariance(const reco::BeamSpot &); + double update(const reco::BeamSpot &, float beam_weight, const bool fill_covariances = false); + void make_vtx_trk_map(const double); + bool clean(); + void remove_vertex(unsigned int); + + // track information + std::vector trackinfo_; + + // vertex lists: + std::vector xv_; + std::vector yv_; + std::vector zv_; + std::vector covv_; + // track-vertex-mapping and weights after a coarse z-cut: + std::vector tkfirstv_; // parallel to the vertex list + std::vector tkmap_; // parallel to tkweight + std::vector tkweight_; // parallel to tkmap + // configuration + double chi_cutoff_; + double z_cutoff_; + double min_trackweight_; + double multivertexfit_; +}; +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h deleted file mode 100644 index d9f7a485f2040..0000000000000 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef DAClusterizerInZ_h -#define DAClusterizerInZ_h - -/**\class DAClusterizerInZ - - Description: separates event tracks into clusters along the beam line - -*/ - -#include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" -#include "TrackingTools/TransientTrack/interface/TransientTrack.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include -#include "DataFormats/Math/interface/Error.h" -#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" -#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" - -class DAClusterizerInZ : public TrackClusterizerInZ { -public: - struct track_t { - double z; // z-coordinate at point of closest approach to the beamline - double dz2; // square of the error of z(pca) - const reco::TransientTrack *tt; // a pointer to the Transient Track - double Z; // Z[i] for DA clustering - double pi; // track weight - }; - - struct vertex_t { - double z; // z coordinate - double pk; // vertex weight for "constrained" clustering - // --- temporary numbers, used during update - double ei; - double sw; - double swz; - double se; - // ---for Tc - double swE; - double Tc; - }; - - DAClusterizerInZ(const edm::ParameterSet &conf); - - std::vector > clusterize( - const std::vector &tracks) const override; - - std::vector vertices(const std::vector &tracks, const int verbosity = 0) const; - - std::vector fill(const std::vector &tracks) const; - - bool split(double beta, std::vector &tks, std::vector &y, double threshold) const; - - double update(double beta, std::vector &tks, std::vector &y) const; - - double update(double beta, std::vector &tks, std::vector &y, double &) const; - - void dump(const double beta, - const std::vector &y, - const std::vector &tks, - const int verbosity = 0) const; - bool merge(std::vector &, int) const; - bool merge(std::vector &, double &) const; - bool purge(std::vector &, std::vector &, double &, const double) const; - - void splitAll(std::vector &y) const; - - double beta0(const double betamax, std::vector &tks, std::vector &y) const; - - double Eik(const track_t &t, const vertex_t &k) const; - -private: - bool verbose_; - bool useTc_; - float vertexSize_; - int maxIterations_; - double coolingFactor_; - float betamax_; - float betastop_; - double dzCutOff_; - double d0CutOff_; -}; - -#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h index b1c61fa9e6f6e..67c3943f9d6d2 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h @@ -11,6 +11,7 @@ #include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" #include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include @@ -264,7 +265,7 @@ class DAClusterizerInZT_vect final : public TrackClusterizerInZ { std::vector > clusterize( const std::vector &tracks) const override; - std::vector vertices(const std::vector &tracks) const; + std::vector vertices(const std::vector &tracks) const override; track_t fill(const std::vector &tracks) const; diff --git a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h index 85fa682a30f85..092ab552e3402 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h +++ b/RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h @@ -175,8 +175,10 @@ class DAClusterizerInZ_vect final : public TrackClusterizerInZ { std::vector > clusterize( const std::vector &tracks) const override; - std::vector vertices(const std::vector &tracks) const; + std::vector vertices(const std::vector &tracks) const override; + std::vector vertices_no_blocks(const std::vector &tracks) const; std::vector vertices_in_blocks(const std::vector &tracks) const; + std::vector fill_vertices(double beta, double rho0, track_t &tracks, vertex_t &vertices) const; track_t fill(const std::vector &tracks) const; diff --git a/RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h b/RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h index d70d2e0291aa7..fa9a127397c2d 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h +++ b/RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h @@ -23,7 +23,7 @@ class GapClusterizerInZ : public TrackClusterizerInZ { float zSeparation() const; - std::vector vertices(const std::vector& tracks) const; + std::vector vertices(const std::vector& tracks) const override; ~GapClusterizerInZ() override{}; diff --git a/RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h b/RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h index 65964e09cfd88..7018aac859c6b 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h +++ b/RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h @@ -39,7 +39,6 @@ class HITrackFilterForPVFinding : public TrackFilterForPVFinding { } static void fillPSetDescription(edm::ParameterSetDescription& desc) { - TrackFilterForPVFinding::fillPSetDescription(desc); desc.add("numTracksThreshold", 0); // HI only desc.add("maxNumTracksThreshold", std::numeric_limits::max()); desc.add("minPtTight", 0.0); diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h new file mode 100644 index 0000000000000..870654e2766df --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h @@ -0,0 +1,33 @@ +#ifndef PrimaryVertexFitterBase_h +#define PrimaryVertexFitterBase_h + +#include + +/**\class PrimaryVertexFitterBase + + Description: base class for primary vertex fitters + +*/ +namespace edm { + class ParameterSet; + class ParameterSetDescription; +} // namespace edm + +namespace reco { + class BeamSpot; + class TransientTrack; +} // namespace reco + +class TransientVertex; + +class PrimaryVertexFitterBase { +public: + PrimaryVertexFitterBase(const edm::ParameterSet &conf) {} + PrimaryVertexFitterBase() {} + virtual ~PrimaryVertexFitterBase() = default; + virtual std::vector fit(const std::vector &, + const std::vector &, + const reco::BeamSpot &, + const bool) = 0; +}; +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h index 8cd9bb51a4ded..fea7068469a59 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h @@ -28,7 +28,6 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" -//#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h" #include "TrackingTools/TransientTrack/interface/TransientTrack.h" #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" #include "TrackingTools/Records/interface/TransientTrackRecord.h" @@ -40,16 +39,23 @@ #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h" -#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h" -#include "RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h" #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" -//#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" +#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h" +#include "RecoVertex/PrimaryVertexProducer/interface/SequentialPrimaryVertexFitterAdapter.h" +#include "RecoVertex/PrimaryVertexProducer/interface/AdaptiveChisquarePrimaryVertexFitter.h" +#include "RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h" + #include "RecoVertex/VertexPrimitives/interface/VertexException.h" #include #include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h" #include "RecoVertex/VertexTools/interface/VertexCompatibleWithBeam.h" #include "DataFormats/Common/interface/ValueMap.h" +// vertex timing +#include "RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmBase.h" +#include "RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h" +#include "RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmLegacy4D.h" + // // class declaration // @@ -75,11 +81,12 @@ class PrimaryVertexProducer : public edm::stream::EDProducer<> { // vtx fitting algorithms struct algo { - VertexFitter<5>* fitter; + PrimaryVertexFitterBase* pv_fitter; VertexCompatibleWithBeam* vertexSelector; std::string label; bool useBeamConstraint; double minNdof; + VertexTimeAlgorithmBase* pv_time_estimator; }; std::vector algorithms; @@ -94,7 +101,11 @@ class PrimaryVertexProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT trkToken; edm::EDGetTokenT > trkTimesToken; edm::EDGetTokenT > trkTimeResosToken; + edm::EDGetTokenT > trackMTDTimeQualityToken; - bool f4D; - bool weightFit; + bool useTransientTrackTime_; + bool useMVASelection_; + edm::ValueMap trackMTDTimeQualities_; + edm::ValueMap trackTimes_; + double minTrackTimeQuality_; }; diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h index e41609be01220..f61844b09aac2 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h @@ -39,10 +39,8 @@ #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h" #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h" -#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h" #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" -//#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" #include "RecoVertex/VertexPrimitives/interface/VertexException.h" #include #include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h" diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexTrackClusterizer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexTrackClusterizer.h new file mode 100644 index 0000000000000..98bddb4bbfd3d --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexTrackClusterizer.h @@ -0,0 +1,28 @@ +#ifndef PrimaryVertexTrackClusterizer_h +#define PrimaryVertexTrackClusterizer_h + +/**\class PrimaryVertexTrackClusterizer + + Description: interface/base class for track clusterizers that separate event tracks into clusters along the beam line + extends TrackClusterizerInZ + +*/ + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" + +class PrimaryVertexTrackClusterizer : public TrackClusterizerInZ { +public: + PrimaryVertexTrackClusterizer() = default; + PrimaryVertexTrackClusterizer(const edm::ParameterSet& conf){}; + virtual std::vector vertices(const std::vector& tracks) const = 0; + virtual std::vector > clusterize( + const std::vector& tracks) const = 0; + + virtual ~PrimaryVertexTrackClusterizer() = default; +}; + +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/SequentialPrimaryVertexFitterAdapter.h b/RecoVertex/PrimaryVertexProducer/interface/SequentialPrimaryVertexFitterAdapter.h new file mode 100644 index 0000000000000..bd5f866e2f21f --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/SequentialPrimaryVertexFitterAdapter.h @@ -0,0 +1,46 @@ +#ifndef SequentialPrimaryVertexFitterAdapter_h +#define SequentialPrimaryVertexFitterAdapter_h + +/**\class SequentialPrimaryVertexFitterAdapter + + Description: Adapter class for Kalman and Adaptive vertex fitters + +*/ + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h" +#include "RecoVertex/VertexPrimitives/interface/VertexFitter.h" + +class SequentialPrimaryVertexFitterAdapter : public PrimaryVertexFitterBase { +public: + SequentialPrimaryVertexFitterAdapter() : fitter(nullptr){}; + SequentialPrimaryVertexFitterAdapter(const VertexFitter<5>* vertex_fitter) : fitter(vertex_fitter){}; + ~SequentialPrimaryVertexFitterAdapter() override = default; + + std::vector fit(const std::vector& dummy, + const std::vector& clusters, + const reco::BeamSpot& beamspot, + const bool useBeamConstraint) override { + std::vector pvs; + for (auto& cluster : clusters) { + const std::vector& tracklist = cluster.originalTracks(); + TransientVertex v; + if (useBeamConstraint && (tracklist.size() > 1)) { + v = fitter->vertex(tracklist, beamspot); + } else if (!(useBeamConstraint) && (tracklist.size() > 1)) { + v = fitter->vertex(tracklist); + } // else: no fit ==> v.isValid()=False + + if (v.isValid()) { + pvs.push_back(v); + } + } + return pvs; + }; + +protected: + // configuration + const VertexFitter<5>* fitter; // Kalman or Adaptive +}; +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h b/RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h index e087aa7215d93..2acf274befb81 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h +++ b/RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h @@ -10,16 +10,17 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include #include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" class TrackClusterizerInZ { public: - TrackClusterizerInZ(){}; + TrackClusterizerInZ() = default; TrackClusterizerInZ(const edm::ParameterSet& conf){}; - + virtual std::vector vertices(const std::vector& tracks) const = 0; virtual std::vector > clusterize( const std::vector& tracks) const = 0; - virtual ~TrackClusterizerInZ(){}; + virtual ~TrackClusterizerInZ() = default; }; #endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmBase.h b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmBase.h new file mode 100644 index 0000000000000..198faa631d0b9 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmBase.h @@ -0,0 +1,56 @@ +#ifndef usercode_PrimaryVertexAnalyzer_VertexTimeAlgorithmBase_h +#define usercode_PrimaryVertexAnalyzer_VertexTimeAlgorithmBase_h +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" + +namespace edm { + class Event; + class EventSetup; + class ParameterSet; + class ParameterSetDescription; + class ConsumesCollector; +} // namespace edm + +class VertexTimeAlgorithmBase { +public: + VertexTimeAlgorithmBase(const edm::ParameterSet& conf, edm::ConsumesCollector& iC) {} + virtual ~VertexTimeAlgorithmBase() = default; + VertexTimeAlgorithmBase(const VertexTimeAlgorithmBase&) = delete; + VertexTimeAlgorithmBase& operator=(const VertexTimeAlgorithmBase&) = delete; + + static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {} + + virtual void setEvent(edm::Event& iEvent, edm::EventSetup const& iSetup) = 0; + + /** + * estimate the vertex time and time uncertainty for transient vertex + * + * returns true when a valid time has been determined, otherwise return false + */ + virtual bool vertexTime(float& vtxTime, float& vtxTimeError, TransientVertex const& vtx) const = 0; + + /** + * replace the vertices in the input vector by new vertices with time coordinates + * determined by the vertexTime method + * this implementation does not alter the weights from the previous fit + * must be overridden to change weights, coordinates, tracklists or to add or remove vertices + */ + virtual void fill_vertex_times(std::vector& pvs) { + for (unsigned int i = 0; i < pvs.size(); i++) { + auto vtx = pvs[i]; + if (vtx.isValid()) { + auto vtxTime(0.f), vtxTimeError(-1.f); + if (vertexTime(vtxTime, vtxTimeError, vtx)) { + auto err = vtx.positionError().matrix4D(); + err(3, 3) = vtxTimeError * vtxTimeError; + auto trkWeightMap3d = vtx.weightMap(); + auto vtx_with_time = TransientVertex( + vtx.position(), vtxTime, err, vtx.originalTracks(), vtx.totalChiSquared(), vtx.degreesOfFreedom()); + vtx_with_time.weightMap(trkWeightMap3d); + pvs[i] = vtx_with_time; + } + } + } + } +}; + +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h new file mode 100644 index 0000000000000..1de2e45e59e6a --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h @@ -0,0 +1,51 @@ + +#ifndef usercode_PrimaryVertexAnalyzer_VertexTimeAlgorithmFromTracksPID_h +#define usercode_PrimaryVertexAnalyzer_VertexTimeAlgorithmFromTracksPID_h + +#include "VertexTimeAlgorithmBase.h" + +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "DataFormats/Common/interface/ValueMap.h" + +class VertexTimeAlgorithmFromTracksPID : public VertexTimeAlgorithmBase { +public: + VertexTimeAlgorithmFromTracksPID(const edm::ParameterSet& conf, edm::ConsumesCollector& iC); + ~VertexTimeAlgorithmFromTracksPID() override = default; + + static void fillPSetDescription(edm::ParameterSetDescription& iDesc); + + void setEvent(edm::Event& iEvent, edm::EventSetup const& iSetup) override; + + bool vertexTime(float& vtxTime, float& vtxTimeError, TransientVertex const& vtx) const override; + +protected: + struct TrackInfo { + double trkWeight; + double trkTimeError; + double trkTimeHyp[3]; + }; + + edm::EDGetTokenT> const trackMTDTimeToken_; + edm::EDGetTokenT> const trackMTDTimeErrorToken_; + edm::EDGetTokenT> const trackMTDTimeQualityToken_; + edm::EDGetTokenT> const trackMTDTofPiToken_; + edm::EDGetTokenT> const trackMTDTofKToken_; + edm::EDGetTokenT> const trackMTDTofPToken_; + + double const minTrackVtxWeight_; + double const minTrackTimeQuality_; + double const probPion_; + double const probKaon_; + double const probProton_; + double const Tstart_; + double const coolingFactor_; + + edm::ValueMap trackMTDTimes_; + edm::ValueMap trackMTDTimeErrors_; + edm::ValueMap trackMTDTimeQualities_; + edm::ValueMap trackMTDTofPi_; + edm::ValueMap trackMTDTofK_; + edm::ValueMap trackMTDTofP_; +}; + +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmLegacy4D.h b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmLegacy4D.h new file mode 100644 index 0000000000000..0d84a2730842d --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmLegacy4D.h @@ -0,0 +1,22 @@ +#ifndef usercode_PrimaryVertexAnalyzer_VertexTimeLegacy4D_h +#define usercode_PrimaryVertexAnalyzer_VertexTimeLegacy4D_h + +#include "VertexTimeAlgorithmBase.h" + +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" + +class VertexTimeAlgorithmLegacy4D : public VertexTimeAlgorithmBase { +public: + VertexTimeAlgorithmLegacy4D(const edm::ParameterSet& conf, edm::ConsumesCollector& iC); + ~VertexTimeAlgorithmLegacy4D() override = default; + + static void fillPSetDescription(edm::ParameterSetDescription& iDesc); + + void setEvent(edm::Event& iEvent, edm::EventSetup const& iSetup) override; + + bool vertexTime(float& vtxTime, float& vtxTimeError, TransientVertex const& vtx) const override; +}; + +#endif diff --git a/RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h b/RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h index 20e4ac03d1e70..881bf2014187c 100644 --- a/RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h +++ b/RecoVertex/PrimaryVertexProducer/interface/WeightedMeanFitter.h @@ -6,6 +6,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/BeamSpot/interface/BeamSpot.h" #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexFitterBase.h" namespace WeightedMeanFitter { @@ -160,7 +161,9 @@ namespace WeightedMeanFitter { dist += std::pow(p.first.z() - z, 2) / (std::pow(wz, 2) + err(2, 2)); chi2 += dist; } - TransientVertex v(GlobalPoint(x, y, z), err, iclus, chi2, (int)ndof_x); + float ndof = + ndof_x > 1 ? (2 * ndof_x - 3) : 0.00001; // ndof_x is actually the number of tracks with non-zero weight + TransientVertex v(GlobalPoint(x, y, z), err, iclus, chi2, ndof); return v; } @@ -455,4 +458,71 @@ namespace WeightedMeanFitter { }; // namespace WeightedMeanFitter +// adapter for the multiprimaryvertexfitter scheme +// this code was originally introduced as part of PrimaryVertexProducer.cc +// by Adriano Di Florio , Giorgio Pizzati et.al. in #39995, then moved here with minor modifications +class WeightedMeanPrimaryVertexEstimator : public PrimaryVertexFitterBase { +public: + WeightedMeanPrimaryVertexEstimator() = default; + ~WeightedMeanPrimaryVertexEstimator() override = default; + + std::vector fit(const std::vector& dummy, + const std::vector& clusters, + const reco::BeamSpot& beamSpot, + const bool useBeamConstraint) override { + std::vector pvs; + std::vector seed(1); + + for (auto& cluster : clusters) { + if (cluster.originalTracks().size() > 1) { + std::vector tracklist = cluster.originalTracks(); + TransientVertex::TransientTrackToFloatMap trkWeightMap; + std::vector> points; + if (useBeamConstraint && (tracklist.size() > 1)) { + for (const auto& itrack : tracklist) { + GlobalPoint p = itrack.stateAtBeamLine().trackStateAtPCA().position(); + GlobalPoint err(itrack.stateAtBeamLine().transverseImpactParameter().error(), + itrack.stateAtBeamLine().transverseImpactParameter().error(), + itrack.track().dzError()); + std::pair p2(p, err); + points.push_back(p2); + } + + TransientVertex v = WeightedMeanFitter::weightedMeanOutlierRejectionBeamSpot(points, tracklist, beamSpot); + if (!v.hasTrackWeight()) { + // if the fitter doesn't provide weights, fill dummy values + TransientVertex::TransientTrackToFloatMap trkWeightMap; + for (const auto& trk : v.originalTracks()) { + trkWeightMap[trk] = 1.; + } + v.weightMap(trkWeightMap); + } + if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError)) + pvs.push_back(v); + } else if (!(useBeamConstraint) && (tracklist.size() > 1)) { + for (const auto& itrack : tracklist) { + GlobalPoint p = itrack.impactPointState().globalPosition(); + GlobalPoint err(itrack.track().dxyError(), itrack.track().dxyError(), itrack.track().dzError()); + std::pair p2(p, err); + points.push_back(p2); + } + + TransientVertex v = WeightedMeanFitter::weightedMeanOutlierRejection(points, tracklist); + if (!v.hasTrackWeight()) { + // if the fitter doesn't provide weights, fill dummy values + TransientVertex::TransientTrackToFloatMap trkWeightMap; + for (const auto& trk : v.originalTracks()) { + trkWeightMap[trk] = 1.; + } + v.weightMap(trkWeightMap); + } + if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError)) + pvs.push_back(v); //FIX with constants + } + } + } + return pvs; + } +}; + #endif diff --git a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc index 6282d76774da7..a196f7e6d2522 100644 --- a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc +++ b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexProducer.cc @@ -1,4 +1,5 @@ #include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/Common/interface/Handle.h" @@ -18,11 +19,11 @@ PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf) : theTTBToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), theConfig(conf) { fVerbose = conf.getUntrackedParameter("verbose", false); + useMVASelection_ = conf.getParameter("useMVACut"); trkToken = consumes(conf.getParameter("TrackLabel")); bsToken = consumes(conf.getParameter("beamSpotLabel")); - f4D = false; - weightFit = false; + useTransientTrackTime_ = false; // select and configure the track selection std::string trackSelectionAlgorithm = @@ -41,53 +42,84 @@ PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf) if (clusteringAlgorithm == "gap") { theTrackClusterizer = new GapClusterizerInZ( conf.getParameter("TkClusParameters").getParameter("TkGapClusParameters")); - } else if (clusteringAlgorithm == "DA") { - theTrackClusterizer = new DAClusterizerInZ( - conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); - } - // provide the vectorized version of the clusterizer, if supported by the build - else if (clusteringAlgorithm == "DA_vect") { + } else if (clusteringAlgorithm == "DA_vect") { theTrackClusterizer = new DAClusterizerInZ_vect( conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); } else if (clusteringAlgorithm == "DA2D_vect") { theTrackClusterizer = new DAClusterizerInZT_vect( conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); - f4D = true; - } - - else { + useTransientTrackTime_ = true; + } else { throw VertexException("PrimaryVertexProducer: unknown clustering algorithm: " + clusteringAlgorithm); } - if (f4D) { - trkTimesToken = consumes>(conf.getParameter("TrackTimesLabel")); - trkTimeResosToken = consumes>(conf.getParameter("TrackTimeResosLabel")); + if (useTransientTrackTime_) { + trkTimesToken = consumes >(conf.getParameter("TrackTimesLabel")); + trkTimeResosToken = consumes >(conf.getParameter("TrackTimeResosLabel")); + trackMTDTimeQualityToken = + consumes >(conf.getParameter("trackMTDTimeQualityVMapTag")); + minTrackTimeQuality_ = conf.getParameter("minTrackTimeQuality"); } // select and configure the vertex fitters std::vector vertexCollections = - conf.getParameter>("vertexCollections"); + conf.getParameter >("vertexCollections"); for (std::vector::const_iterator algoconf = vertexCollections.begin(); algoconf != vertexCollections.end(); algoconf++) { algo algorithm; + + algorithm.label = algoconf->getParameter("label"); + + // configure the fitter and selector std::string fitterAlgorithm = algoconf->getParameter("algorithm"); if (fitterAlgorithm == "KalmanVertexFitter") { - algorithm.fitter = new KalmanVertexFitter(); + algorithm.pv_fitter = new SequentialPrimaryVertexFitterAdapter(new KalmanVertexFitter()); } else if (fitterAlgorithm == "AdaptiveVertexFitter") { - algorithm.fitter = new AdaptiveVertexFitter(GeometricAnnealing(algoconf->getParameter("chi2cutoff"))); + auto fitter = new AdaptiveVertexFitter(GeometricAnnealing(algoconf->getParameter("chi2cutoff"))); + algorithm.pv_fitter = new SequentialPrimaryVertexFitterAdapter(fitter); + } else if (fitterAlgorithm.empty()) { + algorithm.pv_fitter = nullptr; + } else if (fitterAlgorithm == "AdaptiveChisquareVertexFitter") { + algorithm.pv_fitter = new AdaptiveChisquarePrimaryVertexFitter(algoconf->getParameter("chi2cutoff"), + algoconf->getParameter("zcutoff"), + algoconf->getParameter("mintrkweight"), + false); + } else if (fitterAlgorithm == "MultiPrimaryVertexFitter") { + algorithm.pv_fitter = new AdaptiveChisquarePrimaryVertexFitter(algoconf->getParameter("chi2cutoff"), + algoconf->getParameter("zcutoff"), + algoconf->getParameter("mintrkweight"), + true); } else if (fitterAlgorithm == "WeightedMeanFitter") { - algorithm.fitter = nullptr; - weightFit = true; + algorithm.pv_fitter = new WeightedMeanPrimaryVertexEstimator(); } else { throw VertexException("PrimaryVertexProducer: unknown algorithm: " + fitterAlgorithm); } - algorithm.label = algoconf->getParameter("label"); algorithm.minNdof = algoconf->getParameter("minNdof"); algorithm.useBeamConstraint = algoconf->getParameter("useBeamConstraint"); algorithm.vertexSelector = new VertexCompatibleWithBeam(VertexDistanceXY(), algoconf->getParameter("maxDistanceToBeam")); + + // configure separate vertex time reconstruction if applicable + // note that the vertex time could, in principle, also come from the clusterizer or the vertex fit + + const auto& pv_time_conf = algoconf->getParameter("vertexTimeParameters"); + const std::string vertexTimeAlgorithm = pv_time_conf.getParameter("algorithm"); + edm::ConsumesCollector&& collector = consumesCollector(); + + if (vertexTimeAlgorithm.empty()) { + algorithm.pv_time_estimator = nullptr; + } else if (vertexTimeAlgorithm == "legacy4D") { + useTransientTrackTime_ = true; + algorithm.pv_time_estimator = + new VertexTimeAlgorithmLegacy4D(pv_time_conf.getParameter("legacy4D"), collector); + } else if (vertexTimeAlgorithm == "fromTracksPID") { + algorithm.pv_time_estimator = new VertexTimeAlgorithmFromTracksPID( + pv_time_conf.getParameter("fromTracksPID"), collector); + } else { + edm::LogWarning("MisConfiguration") << "unknown vertexTimeParameters.algorithm" << vertexTimeAlgorithm; + } algorithms.push_back(algorithm); produces(algorithm.label); @@ -113,8 +145,10 @@ PrimaryVertexProducer::~PrimaryVertexProducer() { if (theTrackClusterizer) delete theTrackClusterizer; for (std::vector::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) { - if (algorithm->fitter) - delete algorithm->fitter; + if (algorithm->pv_fitter) + delete algorithm->pv_fitter; + if (algorithm->pv_time_estimator) + delete algorithm->pv_time_estimator; if (algorithm->vertexSelector) delete algorithm->vertexSelector; } @@ -135,8 +169,8 @@ void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& i VertexState beamVertexState(beamSpot); if ((beamVertexState.error().cxx() <= 0.) || (beamVertexState.error().cyy() <= 0.) || (beamVertexState.error().czz() <= 0.)) { - validBS = false; edm::LogError("UnusableBeamSpot") << "Beamspot with invalid errors " << beamVertexState.error().matrix(); + validBS = false; } //if this is a recovery iteration, check if we already have a valid PV @@ -160,155 +194,93 @@ void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& i edm::Handle tks; iEvent.getByToken(trkToken, tks); + for (auto& algo : algorithms) { + if (algo.pv_time_estimator) { + algo.pv_time_estimator->setEvent(iEvent, iSetup); + } + } + // interface RECO tracks to vertex reconstruction const auto& theB = &iSetup.getData(theTTBToken); std::vector t_tks; - if (f4D) { - edm::Handle> trackTimesH; - edm::Handle> trackTimeResosH; - iEvent.getByToken(trkTimesToken, trackTimesH); - iEvent.getByToken(trkTimeResosToken, trackTimeResosH); - t_tks = (*theB).build(tks, beamSpot, *(trackTimesH.product()), *(trackTimeResosH.product())); + if (useTransientTrackTime_) { + auto const& trackTimeResos_ = iEvent.get(trkTimeResosToken); + auto trackTimes_ = iEvent.get(trkTimesToken); + + if (useMVASelection_) { + trackMTDTimeQualities_ = iEvent.get(trackMTDTimeQualityToken); + + for (unsigned int i = 0; i < (*tks).size(); i++) { + const reco::TrackRef ref(tks, i); + auto const trkTimeQuality = trackMTDTimeQualities_[ref]; + if (trkTimeQuality < minTrackTimeQuality_) { + trackTimes_[ref] = std::numeric_limits::max(); + } + } + t_tks = (*theB).build(tks, beamSpot, trackTimes_, trackTimeResos_); + } else { + t_tks = (*theB).build(tks, beamSpot, trackTimes_, trackTimeResos_); + } } else { t_tks = (*theB).build(tks, beamSpot); } - if (fVerbose) { - std::cout << "RecoVertex/PrimaryVertexProducer" - << "Found: " << t_tks.size() << " reconstructed tracks" - << "\n"; - } // select tracks std::vector&& seltks = theTrackFilter->select(t_tks); // clusterize tracks in Z - std::vector>&& clusters = theTrackClusterizer->clusterize(seltks); + std::vector&& clusters = theTrackClusterizer->vertices(seltks); if (fVerbose) { - std::cout << " clustering returned " << clusters.size() << " clusters from " << seltks.size() - << " selected tracks" << std::endl; + edm::LogPrint("PrimaryVertexProducer") + << "Clustering returned " << clusters.size() << " clusters from " << seltks.size() << " selected tracks"; } // vertex fits for (std::vector::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) { auto result = std::make_unique(); reco::VertexCollection& vColl = (*result); - std::vector pvs; - for (std::vector>::const_iterator iclus = clusters.begin(); - iclus != clusters.end(); - iclus++) { - double sumwt = 0.; - double sumwt2 = 0.; - double sumw = 0.; - double meantime = 0.; - double vartime = 0.; - if (f4D) { - for (const auto& tk : *iclus) { - const double time = tk.timeExt(); - const double err = tk.dtErrorExt(); - const double inverr = err > 0. ? 1.0 / err : 0.; - const double w = inverr * inverr; - sumwt += w * time; - sumwt2 += w * time * time; - sumw += w; - } - meantime = sumwt / sumw; - double sumsq = sumwt2 - sumwt * sumwt / sumw; - double chisq = iclus->size() > 1 ? sumsq / double(iclus->size() - 1) : sumsq / double(iclus->size()); - vartime = chisq / sumw; - } - TransientVertex v; - if (algorithm->fitter) { - if (algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) { - v = algorithm->fitter->vertex(*iclus, beamSpot); - } else if (!(algorithm->useBeamConstraint) && (iclus->size() > 1)) { - v = algorithm->fitter->vertex(*iclus); - } // else: no fit ==> v.isValid()=False - } else if (weightFit) { - std::vector> points; - if (algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) { - for (const auto& itrack : *iclus) { - GlobalPoint p = itrack.stateAtBeamLine().trackStateAtPCA().position(); - GlobalPoint err(itrack.stateAtBeamLine().transverseImpactParameter().error(), - itrack.stateAtBeamLine().transverseImpactParameter().error(), - itrack.track().dzError()); - std::pair p2(p, err); - points.push_back(p2); - } - - v = WeightedMeanFitter::weightedMeanOutlierRejectionBeamSpot(points, *iclus, beamSpot); - if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError)) - pvs.push_back(v); - } else if (!(algorithm->useBeamConstraint) && (iclus->size() > 1)) { - for (const auto& itrack : *iclus) { - GlobalPoint p = itrack.impactPointState().globalPosition(); - GlobalPoint err(itrack.track().dxyError(), itrack.track().dxyError(), itrack.track().dzError()); - std::pair p2(p, err); - points.push_back(p2); - } - - v = WeightedMeanFitter::weightedMeanOutlierRejection(points, *iclus); - if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError)) - pvs.push_back(v); //FIX with constants - } - } else - throw VertexException( - "PrimaryVertexProducer: Something went wrong. You are not using the weighted mean fit and no algorithm was " - "selected."); - - // 4D vertices: add timing information - if (f4D and v.isValid()) { - auto err = v.positionError().matrix4D(); - err(3, 3) = vartime; - auto trkWeightMap3d = v.weightMap(); // copy the 3d-fit weights - v = TransientVertex(v.position(), meantime, err, v.originalTracks(), v.totalChiSquared(), v.degreesOfFreedom()); - v.weightMap(trkWeightMap3d); - } + if (algorithm->pv_fitter == nullptr) { + pvs = clusters; + } else { + pvs = algorithm->pv_fitter->fit(seltks, clusters, beamSpot, algorithm->useBeamConstraint); + } + + if (algorithm->pv_time_estimator != nullptr) { + algorithm->pv_time_estimator->fill_vertex_times(pvs); + } + + // sort vertices by pt**2 vertex + if (pvs.size() > 1) { + sort(pvs.begin(), pvs.end(), VertexHigherPtSquared()); + } - if (fVerbose) { - if (v.isValid()) { - std::cout << "x,y,z"; - if (f4D) - std::cout << ",t"; - std::cout << "=" << v.position().x() << " " << v.position().y() << " " << v.position().z(); - if (f4D) - std::cout << " " << v.time(); - std::cout << " cluster size = " << (*iclus).size() << std::endl; - } else { - std::cout << "Invalid fitted vertex, cluster size=" << (*iclus).size() << std::endl; + // select and convert transient vertices to (reco) vertices + for (std::vector::const_iterator iv = pvs.begin(); iv != pvs.end(); iv++) { + if (iv->isValid() && (iv->degreesOfFreedom() >= algorithm->minNdof)) { + reco::Vertex v = *iv; + if (!validBS || ((*(algorithm->vertexSelector))(v, beamVertexState))) { + vColl.push_back(v); } } - - //for weightFit we have already pushed it above (no timing infomration anyway) - if (v.isValid() && not weightFit && (v.degreesOfFreedom() >= algorithm->minNdof) && - (!validBS || (*(algorithm->vertexSelector))(v, beamVertexState))) - pvs.push_back(v); - } // end of cluster loop + } if (fVerbose) { - std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << pvs.size() << std::endl; + edm::LogPrint("PrimaryVertexProducer") << "PrimaryVertexProducer \"" << algorithm->label << "\" contains " + << pvs.size() << " reco::Vertex candidates"; } - if (clusters.size() > 2 && clusters.size() > 2 * pvs.size()) + if (clusters.size() > 2 && clusters.size() > 2 * pvs.size()) { edm::LogWarning("PrimaryVertexProducer") - << "more than half of candidate vertices lost " << pvs.size() << ' ' << clusters.size(); - - if (pvs.empty() && seltks.size() > 5) - edm::LogWarning("PrimaryVertexProducer") - << "no vertex found with " << seltks.size() << " tracks and " << clusters.size() << " vertex-candidates"; - - // sort vertices by pt**2 vertex (aka signal vertex tagging) - if (pvs.size() > 1) { - sort(pvs.begin(), pvs.end(), VertexHigherPtSquared()); + << "More than 50% of candidate vertices lost (" << pvs.size() << " out of " << clusters.size() << ")"; } - // convert transient vertices returned by the theAlgo to (reco) vertices - for (std::vector::const_iterator iv = pvs.begin(); iv != pvs.end(); iv++) { - reco::Vertex v = *iv; - vColl.push_back(v); + if (pvs.empty() && seltks.size() > 5) { + edm::LogWarning("PrimaryVertexProducer") + << "No vertex found with " << seltks.size() << " tracks and " << clusters.size() << " vertex candidates"; } if (vColl.empty()) { @@ -319,16 +291,14 @@ void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& i we(1, 1) = 10000; we(2, 2) = 10000; vColl.push_back(reco::Vertex(beamSpot.position(), we, 0., 0., 0)); - if (fVerbose) { - std::cout << "RecoVertex/PrimaryVertexProducer: " - << "Beamspot with invalid errors " << bse.matrix() << std::endl; - std::cout << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n"; - } + edm::LogWarning("PrimaryVertexProducer") << "Zero recostructed vertices, will put reco::Vertex derived from " + "dummy/fake BeamSpot into Event, BeamSpot has invalid errors: " + << bse.matrix(); } else { vColl.push_back(reco::Vertex(beamSpot.position(), beamSpot.rotatedCovariance3D(), 0., 0., 0)); if (fVerbose) { - std::cout << "RecoVertex/PrimaryVertexProducer: " - << " will put Vertex derived from BeamSpot into Event.\n"; + edm::LogWarning("PrimaryVertexProducer") + << "Zero recostructed vertices, will put reco::Vertex derived from BeamSpot into Event."; } } } @@ -336,15 +306,18 @@ void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& i if (fVerbose) { int ivtx = 0; for (reco::VertexCollection::const_iterator v = vColl.begin(); v != vColl.end(); ++v) { - std::cout << "recvtx " << ivtx++ << "#trk " << std::setw(3) << v->tracksSize() << " chi2 " << std::setw(4) - << v->chi2() << " ndof " << std::setw(3) << v->ndof() << " x " << std::setw(6) << v->position().x() - << " dx " << std::setw(6) << v->xError() << " y " << std::setw(6) << v->position().y() << " dy " - << std::setw(6) << v->yError() << " z " << std::setw(6) << v->position().z() << " dz " << std::setw(6) - << v->zError(); - if (f4D) { - std::cout << " t " << std::setw(6) << v->t() << " dt " << std::setw(6) << v->tError(); + edm::LogPrint("PrimaryVertexProducer") + << "recvtx " << std::setw(3) << std::fixed << ivtx++ << " #trk " << std::setw(3) << v->tracksSize() + << " chi2 " << std::setw(5) << std::setprecision(1) << v->chi2() << " ndof " << std::setw(5) + << std::setprecision(1) << v->ndof() << " x " << std::setw(7) << std::setprecision(4) << v->position().x() + << " dx " << std::setw(6) << std::setprecision(4) << v->xError() << " y " << std::setw(7) + << std::setprecision(4) << v->position().y() << " dy " << std::setw(6) << std::setprecision(4) + << v->yError() << " z " << std::setw(8) << std::setprecision(4) << v->position().z() << " dz " + << std::setw(6) << std::setprecision(4) << v->zError(); + if (v->tError() > 0) { + edm::LogPrint("PrimaryVertexProducer") << " t " << std::setw(6) << std::setprecision(3) << v->t() << " dt " + << std::setw(6) << std::setprecision(3) << v->tError(); } - std::cout << std::endl; } } @@ -353,7 +326,19 @@ void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& i } void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - // offlinePrimaryVertices + edm::ParameterSetDescription psd_pv_time; + { + edm::ParameterSetDescription psd1; + VertexTimeAlgorithmFromTracksPID::fillPSetDescription(psd1); + psd_pv_time.add("fromTracksPID", psd1); + + edm::ParameterSetDescription psd2; + VertexTimeAlgorithmLegacy4D::fillPSetDescription(psd2); + psd_pv_time.add("legacy4D", psd2); + } + psd_pv_time.add("algorithm", ""); // default = none + + // vertex collections edm::ParameterSetDescription desc; { edm::ParameterSetDescription vpsd1; @@ -362,7 +347,12 @@ void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& des vpsd1.add("useBeamConstraint", false); vpsd1.add("label", ""); vpsd1.add("chi2cutoff", 2.5); + vpsd1.add("zcutoff", 1.0); + vpsd1.add("mintrkweight", 0.0); vpsd1.add("minNdof", 0.0); + vpsd1.add("vertexTimeParameters", psd_pv_time); + + // two default values : with- and without beam constraint std::vector temp1; temp1.reserve(2); { @@ -372,7 +362,12 @@ void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& des temp2.addParameter("useBeamConstraint", false); temp2.addParameter("label", ""); temp2.addParameter("chi2cutoff", 2.5); + temp2.addParameter("zcutoff", 1.0); + temp2.addParameter("mintrkweight", 0.); temp2.addParameter("minNdof", 0.0); + edm::ParameterSet temp_vertexTime; + temp_vertexTime.addParameter("algorithm", ""); + temp2.addParameter("vertexTimeParameters", temp_vertexTime); temp1.push_back(temp2); } { @@ -382,7 +377,12 @@ void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& des temp2.addParameter("useBeamConstraint", true); temp2.addParameter("label", "WithBS"); temp2.addParameter("chi2cutoff", 2.5); + temp2.addParameter("zcutoff", 1.0); + temp2.addParameter("mintrkweight", 0.); temp2.addParameter("minNdof", 2.0); + edm::ParameterSet temp_vertexTime; + temp_vertexTime.addParameter("algorithm", ""); + temp2.addParameter("vertexTimeParameters", temp_vertexTime); temp1.push_back(temp2); } desc.addVPSet("vertexCollections", vpsd1, temp1); @@ -391,15 +391,14 @@ void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& des { edm::ParameterSetDescription psd0; TrackFilterForPVFinding::fillPSetDescription(psd0); - psd0.add("numTracksThreshold", 0); // HI only - psd0.add("maxNumTracksThreshold", 10000000); // HI only - psd0.add("minPtTight", 0.0); // HI only + HITrackFilterForPVFinding::fillPSetDescription(psd0); // HI only desc.add("TkFilterParameters", psd0); } desc.add("beamSpotLabel", edm::InputTag("offlineBeamSpot")); desc.add("TrackLabel", edm::InputTag("generalTracks")); - desc.add("TrackTimeResosLabel", edm::InputTag("dummy_default")); // 4D only - desc.add("TrackTimesLabel", edm::InputTag("dummy_default")); // 4D only + desc.add("TrackTimeResosLabel", edm::InputTag("dummy_default")); // 4D only + desc.add("TrackTimesLabel", edm::InputTag("dummy_default")); // 4D only + desc.add("trackMTDTimeQualityVMapTag", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA")); // 4D only { edm::ParameterSetDescription psd0; @@ -418,6 +417,8 @@ void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& des desc.add("isRecoveryIteration", false); desc.add("recoveryVtxCollection", {""}); + desc.add("useMVACut", false); + desc.add("minTrackTimeQuality", 0.8); descriptions.add("primaryVertexProducer", desc); } diff --git a/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVertices_cfi.py b/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVertices_cfi.py index ec9a05c3a4393..cff8dbaebaa1e 100644 --- a/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVertices_cfi.py +++ b/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVertices_cfi.py @@ -1,51 +1,45 @@ import FWCore.ParameterSet.Config as cms -from RecoVertex.PrimaryVertexProducer.TkClusParameters_cff import DA_vectParameters +from RecoVertex.PrimaryVertexProducer.primaryVertexProducer_cfi import primaryVertexProducer -offlinePrimaryVertices = cms.EDProducer( - "PrimaryVertexProducer", +offlinePrimaryVertices = primaryVertexProducer.clone() - verbose = cms.untracked.bool(False), - TrackLabel = cms.InputTag("generalTracks"), - beamSpotLabel = cms.InputTag("offlineBeamSpot"), - - TkFilterParameters = cms.PSet( - algorithm=cms.string('filter'), - maxNormalizedChi2 = cms.double(10.0), - minPixelLayersWithHits=cms.int32(2), - minSiliconLayersWithHits = cms.int32(5), - maxD0Significance = cms.double(4.0), - maxD0Error = cms.double(1.0), - maxDzError = cms.double(1.0), - minPt = cms.double(0.0), - maxEta = cms.double(2.4), - trackQuality = cms.string("any") - ), +DA_vectParameters = cms.PSet(primaryVertexProducer.TkClusParameters.clone()) - TkClusParameters = DA_vectParameters, +from Configuration.ProcessModifiers.vertexInBlocks_cff import vertexInBlocks +vertexInBlocks.toModify(offlinePrimaryVertices, + TkClusParameters = dict( + TkDAClusParameters = dict( + runInBlocks = True, + block_size = 128, + overlap_frac = 0.5 + ) + ) +) - vertexCollections = cms.VPSet( - [cms.PSet(label=cms.string(""), - algorithm=cms.string("AdaptiveVertexFitter"), - chi2cutoff = cms.double(2.5), - minNdof=cms.double(0.0), - useBeamConstraint = cms.bool(False), - maxDistanceToBeam = cms.double(1.0) - ), - cms.PSet(label=cms.string("WithBS"), - algorithm = cms.string('AdaptiveVertexFitter'), - chi2cutoff = cms.double(2.5), - minNdof=cms.double(2.0), - useBeamConstraint = cms.bool(True), - maxDistanceToBeam = cms.double(1.0), - ) - ] - ), - - isRecoveryIteration = cms.bool(False), - recoveryVtxCollection = cms.InputTag("") +from Configuration.Eras.Modifier_phase2_tracker_cff import phase2_tracker +(phase2_tracker & vertexInBlocks).toModify(offlinePrimaryVertices, + TkClusParameters = dict( + TkDAClusParameters = dict( + block_size = 512, + overlap_frac = 0.5) + ) +) - +from Configuration.Eras.Modifier_highBetaStar_2018_cff import highBetaStar_2018 +highBetaStar_2018.toModify(offlinePrimaryVertices, + TkClusParameters = dict( + TkDAClusParameters = dict( + Tmin = 4.0, + Tpurge = 1.0, + Tstop = 1.0, + vertexSize = 0.01, + d0CutOff = 4., + dzCutOff = 5., + zmerge = 2.e-2, + uniquetrkweight = 0.9 + ) + ) ) from Configuration.ProcessModifiers.weightedVertexing_cff import weightedVertexing @@ -97,11 +91,8 @@ maxNumTracksThreshold = cms.int32(1000), minPtTight = cms.double(1.0) ), - TkClusParameters = cms.PSet( - algorithm = cms.string("gap"), - TkGapClusParameters = cms.PSet( - zSeparation = cms.double(1.0) - ) + TkClusParameters = dict( + algorithm = "gap" ) ) @@ -121,4 +112,3 @@ 1: dict(chi2cutoff = 4.0, minNdof = -2.0), } ) - diff --git a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py b/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py deleted file mode 100644 index a9fd50d4301dd..0000000000000 --- a/RecoVertex/PrimaryVertexProducer/python/TkClusParameters_cff.py +++ /dev/null @@ -1,77 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -DA_vectParameters = cms.PSet( - algorithm = cms.string("DA_vect"), - TkDAClusParameters = cms.PSet( - coolingFactor = cms.double(0.6), # moderate annealing speed - zrange = cms.double(4.), # consider only clusters within 4 sigma*sqrt(T) of a track - delta_highT = cms.double(1.e-2), # convergence requirement at high T - delta_lowT = cms.double(1.e-3), # convergence requirement at low T - convergence_mode = cms.int32(0), # 0 = two steps, 1 = dynamic with sqrt(T) - Tmin = cms.double(2.0), # end of vertex splitting - Tpurge = cms.double(2.0), # cleaning - Tstop = cms.double(0.5), # end of annealing - vertexSize = cms.double(0.006), # added in quadrature to track-z resolutions - d0CutOff = cms.double(3.), # downweight high IP tracks - dzCutOff = cms.double(3.), # outlier rejection after freeze-out (T= 0) && " negative chi**2"); +#endif + return chsq; +} + +void AdaptiveChisquarePrimaryVertexFitter::fill_trackinfo(const std::vector &tracks, + const reco::BeamSpot &beamSpot) { + /* fill track information used during fits into arrays, parallell to the list of input tracks */ + + trackinfo_.clear(); + trackinfo_.reserve(tracks.size()); + + for (auto &trk : tracks) { + TrackInfo ti; + // F1,F2 are the perigee parameters (3,4) + const auto tspca = trk.stateAtBeamLine().trackStateAtPCA(); // freeTrajectoryState + const auto tspca_pe = PerigeeConversions::ftsToPerigeeError(tspca); + const auto momentum = tspca.momentum(); + auto const cos_phi = momentum.x() / momentum.perp(); + auto const sin_phi = momentum.y() / momentum.perp(); + auto const tan_lambda = momentum.z() / momentum.perp(); + + // covariance matrix of (F1,F2) + double cov11 = tspca_pe.covarianceMatrix()(3, 3); + double cov22 = tspca_pe.covarianceMatrix()(4, 4); + double cov12 = tspca_pe.covarianceMatrix()(3, 4); + + // S = cov^{-1} + double DetV = cov11 * cov22 - cov12 * cov12; + if (fabs(DetV) < 1.e-16) { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") + << "Warning, det(V) almost vanishes : " << DetV << " !! This should not happen!" << std::endl; + ti.S11 = 0; + ti.S22 = 0; + ti.S12 = 0; + } else { + ti.S11 = cov22 / DetV; + ti.S22 = cov11 / DetV; + ti.S12 = -cov12 / DetV; + } + ti.b1 = tspca.position().x() * sin_phi - tspca.position().y() * cos_phi; + ti.H1[0] = -sin_phi; + ti.H1[1] = cos_phi; + ti.H1[2] = 0; + ti.b2 = tspca.position().z() - (tspca.position().x() * cos_phi + tspca.position().y() * sin_phi) * tan_lambda; + ti.H2[0] = cos_phi * tan_lambda; + ti.H2[1] = sin_phi * tan_lambda; + ti.H2[2] = -1.; + + for (int k = 0; k < 3; k++) { + double SH1k = (ti.S11 * ti.H1[k] + ti.S12 * ti.H2[k]); + double SH2k = (ti.S12 * ti.H1[k] + ti.S22 * ti.H2[k]); + ti.g[k] = ti.b1 * SH1k + ti.b2 * SH2k; + for (int l = 0; l < 3; l++) { + ti.C(l, k) = ti.H1[l] * SH1k + ti.H2[l] * SH2k; + } + } + + ti.zpca = tspca.position().z(); + ti.dzError = trk.track().dzError(); + trackinfo_.push_back(ti); + } +} + +void AdaptiveChisquarePrimaryVertexFitter::make_vtx_trk_map(double zrange_scale) { + unsigned const int nv = xv_.size(); + unsigned const int nt = trackinfo_.size(); + +#ifdef PVTX_DEBUG + if (nv < 1) { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFit") << " empty vertex list with " << nt << " tracks" << std::endl; + return; + } +#endif + + // parallel lists for track to vertex mapping, tracks are not sorted + tkmap_.clear(); // index in trackinfo_ + tkweight_.clear(); // weight in vertex + tkfirstv_.clear(); // each vertex k owns a section of those list : tkfirstv_[k] .. tkfirstv_[k+1]-1 + + if (nv == 1) { + // always accept all tracks for a single vertex fit + tkfirstv_.push_back(0); + tkfirstv_.push_back(nt); + tkweight_.assign(nt, 0.); + tkmap_.reserve(nt); + for (unsigned int i = 0; i < nt; i++) { + tkmap_.emplace_back(i); + } + return; + } + + // n > 1 + tkmap_.reserve(nv * 100); + tkweight_.reserve(nv * 100); + for (unsigned int k = 0; k < nv; k++) { + tkfirstv_.emplace_back(tkmap_.size()); + for (unsigned int i = 0; i < nt; i++) { + auto &ti = trackinfo_[i]; + const double zrange = zrange_scale * ti.dzError; + if (std::abs(zv_[k] - ti.zpca) < z_cutoff_) { + const double dztrk = ti.b2 + xv_[k] * ti.H2[0] + yv_[k] * ti.H2[1] - zv_[k]; + if (std::abs(dztrk) < zrange) { + tkmap_.emplace_back(i); + tkweight_.emplace_back(0.); + } + } + } + } + tkfirstv_.emplace_back(tkmap_.size()); // extra entry, simplifies loops, every vertex has a "successor" now +} + +void AdaptiveChisquarePrimaryVertexFitter::fill_weights(const reco::BeamSpot &beamspot, double beta) { + // multi-vertex version + unsigned const int nt = trackinfo_.size(); + unsigned const int nv = xv_.size(); + const double beta_over_2 = 0.5 * beta; + const double argmax = beta_over_2 * chi_cutoff_ * chi_cutoff_ * 5; + const double Z_cutoff = vdt::fast_exp(-beta_over_2 * chi_cutoff_ * chi_cutoff_); + + std::vector Z_track(nt, Z_cutoff); + + // evaluate and cache track-vertex assignment chi**2 for all clusters and sum up Z + for (unsigned int k = 0; k < nv; k++) { + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) { + const unsigned int i = tkmap_[j]; + double arg = beta_over_2 * track_in_vertex_chsq(trackinfo_[i], xv_[k], yv_[k], zv_[k]); + if (arg < argmax) { + const double e = vdt::fast_exp(-arg); + tkweight_[j] = e; // must later be normalized by the proper Z_track[i] + Z_track[i] += e; // sum up exponentials to normalize + } else { + tkweight_[j] = 0.; + } + } + } + + // now we have the partition function, Z_i and can evaluate assignment probabilities (aka weights) + for (unsigned int j = 0; j < tkmap_.size(); j++) { + const unsigned int i = tkmap_[j]; +#ifdef PVT_DEBUG + assert((i < nt) && "tkmap out of range"); + assert((tkmap_.size() == tkweight_.size()) && "map and list not aliged"); +#endif + tkweight_[j] /= Z_track[i]; + } +} + +bool AdaptiveChisquarePrimaryVertexFitter::clean() { + /* in multi-vertex fitting, nearby vertices can fall on top of each other, + even when the initial seeds don't, some kind of duplicate removal is required + the approach in this method is similar to the method applied in clustering: + at least two tracks with a weight above a threshold (trkweight_threshold) are required. + vertices that don't fulfill this are either insignficant or very close + to another vertex + */ + const double trkweight_threshold = 0.7; + unsigned int nv = xv_.size(); + if (nv < 2) + return false; + + // sum of weights per vertex + std::vector wsumhi(nv, 0); + for (unsigned int k = 0; k < nv; k++) { + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) { + if (tkweight_[j] > trkweight_threshold) + wsumhi[k] += tkweight_[j]; + } + } + + double dzmin = 0; + unsigned int k_dzmin = 0; + for (unsigned int k = 0; k < nv - 1; k++) { + double dz = std::abs(zv_[k + 1] - zv_[k]); + if ((k == 0) || (dz < dzmin)) { + dzmin = dz; + k_dzmin = k; + } + } + + if ((std::abs(dzmin) < 0.0200) && (std::min(wsumhi[k_dzmin], wsumhi[k_dzmin + 1]) < 0.5)) { + if (wsumhi[k_dzmin] < wsumhi[k_dzmin + 1]) { + remove_vertex(k_dzmin); + } else { + remove_vertex(k_dzmin + 1); + } + } + + return true; +} + +void AdaptiveChisquarePrimaryVertexFitter::remove_vertex(unsigned int k) { + // remove a vertex or rather merge it with it's neighbour + // used for multi-vertex fits only + unsigned int nv = xv_.size(); + if (nv < 2) + return; + + // 1) remove the vertex from the vertex list + xv_.erase(xv_.begin() + k); + yv_.erase(yv_.begin() + k); + zv_.erase(zv_.begin() + k); + covv_.erase(covv_.begin() + k); + + // 2) adjust the track-map map + // 2a) remove the map entries that belong the deleted vertex + const unsigned int num_erased_map_entries = tkfirstv_[k + 1] - tkfirstv_[k]; + tkmap_.erase(tkmap_.begin() + tkfirstv_[k], tkmap_.begin() + tkfirstv_[k + 1]); + tkweight_.erase(tkweight_.begin() + tkfirstv_[k], tkweight_.begin() + tkfirstv_[k + 1]); + // 2b) adjust pointers for the following vertices, including the dummy entry behind the last (now [nv-1]) + for (unsigned int k1 = k + 1; k1 < nv + 1; k1++) { + tkfirstv_[k1] -= num_erased_map_entries; + } + // 2c) erase the pointer of the removed vertex + tkfirstv_.erase(tkfirstv_.begin() + k); +} + +double AdaptiveChisquarePrimaryVertexFitter::update(const reco::BeamSpot &beamspot, + const float beam_weight, + const bool fill_covariances) { + double rho_vtx = 0; + double delta_z = 0; + double delta_x = 0; + double delta_y = 0; + unsigned const int nt = trackinfo_.size(); + unsigned const int nv = xv_.size(); + if (fill_covariances) { + covv_.clear(); + } + + // initial value for S, 0 or inverse of the beamspot covariance matrix + Error3 S0; + double c_beam[3] = {0, 0, 0}; + if (beam_weight > 0) { + S0 = get_inverse_beam_covariance(beamspot); + for (unsigned int j = 0; j < 3; j++) { + c_beam[j] = -(S0(j, 0) * beamspot.x0() + S0(j, 1) * beamspot.y0() + S0(j, 2) * beamspot.z0()); + } + } + + for (unsigned int k = 0; k < nv; k++) { + rho_vtx = 0; + Error3 S(S0); + // sum track contributions + double c[3] = {c_beam[0], c_beam[1], c_beam[2]}; + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) { + const unsigned int i = tkmap_[j]; + const auto w = tkweight_[j]; + rho_vtx += w; + S += w * trackinfo_[i].C; + for (unsigned int l = 0; l < 3; l++) { + c[l] += w * trackinfo_[i].g[l]; + } + } + +#ifdef PVTX_DEBUG + if ((fabs(S(1, 2) - S(2, 1)) > 1e-3) || (fabs(S(0, 2) - S(2, 0)) > 1e-3) || (fabs(S(0, 1) - S(1, 0)) > 1e-3) || + (S(0, 0) <= 0) || (S(0, 0) <= 0) || (S(0, 0) <= 0)) { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") << "update() bad S-matrix S=" << std::endl + << S << std::endl; + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") + << "n-vertex = " << nv << " n-track = " << nt << std::endl; + } +#endif + + const auto xold = xv_[k]; + const auto yold = yv_[k]; + const auto zold = zv_[k]; + + if (S.Invert()) { + xv_[k] = -(S(0, 0) * c[0] + S(0, 1) * c[1] + S(0, 2) * c[2]); + yv_[k] = -(S(1, 0) * c[0] + S(1, 1) * c[1] + S(1, 2) * c[2]); + zv_[k] = -(S(2, 0) * c[0] + S(2, 1) * c[1] + S(2, 2) * c[2]); + if (fill_covariances) { + covv_.emplace_back(S); + } + } else { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") << "update() Matrix inversion failed" << S << std::endl; + if (fill_covariances) { + Error3 covv_dummy; + covv_dummy(0, 0) = 100.; + covv_dummy(1, 1) = 100.; + covv_dummy(2, 2) = 100.; + covv_.emplace_back(covv_dummy); + } + } + + if ((nt > 1) && (rho_vtx > 1.0)) { + delta_x = std::max(delta_x, std::abs(xv_[k] - xold)); + delta_y = std::max(delta_y, std::abs(yv_[k] - yold)); + delta_z = std::max(delta_z, std::abs(zv_[k] - zold)); + } + + } // vertex loop + + return std::max(delta_z, std::max(delta_x, delta_y)); +} + +AdaptiveChisquarePrimaryVertexFitter::Error3 AdaptiveChisquarePrimaryVertexFitter::get_inverse_beam_covariance( + const reco::BeamSpot &beamspot) { + auto SBeam = beamspot.rotatedCovariance3D(); + if (SBeam.Invert()) { + return SBeam; + } else { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") + << "Warning, beam-spot covariance matrix inversion failed " << std::endl; + Error3 S0; + S0(0, 0) = 1. / pow(beamspot.BeamWidthX(), 2); + S0(1, 1) = 1. / pow(beamspot.BeamWidthY(), 2); + S0(2, 2) = 1. / pow(beamspot.sigmaZ(), 2); + return S0; + } +} + +TransientVertex AdaptiveChisquarePrimaryVertexFitter::get_TransientVertex( + const unsigned int k, + const std::vector> &vertex_track_weights, + const std::vector &tracks, + const float beam_weight, + const reco::BeamSpot &beamspot) { + const GlobalPoint pos(xv_[k], yv_[k], zv_[k]); + const GlobalError posError( + covv_[k](0, 0), covv_[k](1, 0), covv_[k](1, 1), covv_[k](2, 0), covv_[k](2, 1), covv_[k](2, 2)); + float chi2 = 0.; + float vtx_ndof = -3.; + if (beam_weight > 0) { + // add beam-spot chi**2 and degrees of freedom + vtx_ndof = 3 * beam_weight; + const auto S = get_inverse_beam_covariance(beamspot); + const double dx = xv_[k] - beamspot.x0(); + const double dy = yv_[k] - beamspot.y0(); + const double dz = zv_[k] - beamspot.z0(); + chi2 = beam_weight * (S(0, 0) * dx * dx + S(1, 1) * dy * dy + 2 * S(0, 1) * dx * dy + S(2, 2) * dz * dz + + 2 * S(0, 2) * dx * dz + 2 * S(1, 2) * dy * dz); + } + + std::vector vertex_tracks; + TransientVertex::TransientTrackToFloatMap trkWeightMap; + for (const auto &tk : vertex_track_weights) { + const unsigned int i = tk.first; + const float track_weight = tk.second; + if (track_weight >= min_trackweight_) { + vertex_tracks.emplace_back(tracks[i]); + trkWeightMap[tracks[i]] = track_weight; + vtx_ndof += 2 * track_weight; + chi2 += track_weight * track_in_vertex_chsq(trackinfo_[i], xv_[k], yv_[k], zv_[k]); + } + } + + auto vtx = TransientVertex(pos, posError, vertex_tracks, chi2, vtx_ndof); + vtx.weightMap(trkWeightMap); + + return vtx; +} + +std::vector AdaptiveChisquarePrimaryVertexFitter::vertices( + const std::vector &tracks, + const std::vector &clusters, + const reco::BeamSpot &beamspot, + const bool useBeamConstraint) { + // simultaneous fit of all vertices in the input list + + const int max_iterations = 50; + + // initialize the vertices + const unsigned int nv = clusters.size(); + xv_.clear(); + xv_.reserve(nv); + yv_.clear(); + yv_.reserve(nv); + zv_.clear(); + zv_.reserve(nv); + tkfirstv_.clear(); + tkfirstv_.reserve(nv + 1); + covv_.clear(); + covv_.reserve(nv); + + // seeds + for (auto &clu : clusters) { + const double zclu = clu.position().z(); + xv_.emplace_back(beamspot.x(zclu)); + yv_.emplace_back(beamspot.y(zclu)); + zv_.emplace_back(zclu); + } + + fill_trackinfo(tracks, beamspot); + + make_vtx_trk_map(5.); // use tracks within 5 sigma windows (if that is less than z_cutoff_) + + float beam_weight = useBeamConstraint ? 1. : 0.; + + double delta = 0; + unsigned int nit = 0; + while ((nit == 0) || ((delta > 0.0001) && (nit < max_iterations))) { + fill_weights(beamspot); + delta = update(beamspot, beam_weight, false); + nit++; + } + if ((nit >= max_iterations) && (delta > 0.01)) { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") + << "iteration limit reached " << nit << " last delta = " << delta << std::endl + << " nv = " << nv << " nt = " << tracks.size() << std::endl; + } + + // may need to remove collapsed vertices + nit = 0; + while ((xv_.size() > 1) && (nit < max_iterations) && (clean())) { + fill_weights(beamspot); + update(beamspot, beam_weight, false); + nit++; + } + + // fill the covariance matrices + update(beamspot, beam_weight, true); + + // assign tracks to vertices + std::vector track_to_vertex(trackinfo_.size(), nv); + // for each track identify the vertex that wants it most + std::vector maxweight(trackinfo_.size(), -1.); + for (unsigned int k = 0; k < nv; k++) { + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) { + const unsigned int i = tkmap_[j]; + if (tkweight_[j] > maxweight[i]) { + maxweight[i] = tkweight_[j]; + track_to_vertex[i] = k; + } + } + } + + // fill the fit result into transient vertices + std::vector pvs; + for (unsigned int k = 0; k < xv_.size(); k++) { + std::vector> vertex_tracks_weights; + for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) { + unsigned int i = tkmap_[j]; + if (track_to_vertex[i] == k) { + vertex_tracks_weights.emplace_back(tkmap_[j], tkweight_[j]); + } + } + pvs.emplace_back(get_TransientVertex(k, vertex_tracks_weights, tracks, beam_weight, beamspot)); + } + + return pvs; +} + +TransientVertex AdaptiveChisquarePrimaryVertexFitter::refit(const TransientVertex &cluster, + const reco::BeamSpot &beamspot, + const bool useBeamConstraint) { + // fit a single vertex using all tracks in the tracklist + const unsigned int nt = cluster.originalTracks().size(); + const int max_iterations = 50; + + // initialize, vectors with size=1 here to avoid code duplication from the multivertex case in update() + const double zclu = cluster.position().z(); + xv_ = {beamspot.x(zclu)}; + yv_ = {beamspot.y(zclu)}; + zv_ = {zclu}; + tkfirstv_ = {0, nt}; + covv_.clear(); + + fill_trackinfo(cluster.originalTracks(), beamspot); + tkweight_.assign(nt, 0.); + tkmap_.clear(); + tkmap_.reserve(nt); + for (unsigned int i = 0; i < nt; i++) { + tkmap_.emplace_back(i); // trivial map for single vertex fits + } + + float beam_weight = useBeamConstraint ? 1. : 0.; + + double delta = 0; + unsigned int nit = 0; + while ((nit == 0) || ((delta > 0.0001) && (nit < max_iterations))) { + fill_weights(beamspot); + delta = update(beamspot, beam_weight, false); + nit++; + } + + if ((nit >= max_iterations) && (delta > 0.1)) { + edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") + << "single vertex fit, iteration limit reached " << nit << " last delta = " << delta << std::endl + << " nt = " << cluster.originalTracks().size() << std::endl; + } + + // fill the covariance matrices + update(beamspot, beam_weight, true); + + // put the result into a transient vertex + std::vector> vertex_track_weights; + for (unsigned int i = 0; i < nt; i++) { + vertex_track_weights.emplace_back(i, tkweight_[i]); + } + + return get_TransientVertex(0, vertex_track_weights, cluster.originalTracks(), beam_weight, beamspot); +} + +// +std::vector AdaptiveChisquarePrimaryVertexFitter::fit(const std::vector &tracks, + const std::vector &clusters, + const reco::BeamSpot &beamspot, + const bool useBeamConstraint) { + if (multivertexfit_) { + return vertices(tracks, clusters, beamspot, useBeamConstraint); + + } else { + // fit the clusters one-by-one using the tracklist of the clusters (ignores the "tracks" argument) + std::vector pvs; + pvs.reserve(clusters.size()); + for (auto &cluster : clusters) { + if (cluster.originalTracks().size() > (useBeamConstraint ? 0 : 1)) { + auto pv = refit(cluster, beamspot, useBeamConstraint); + if (pv.isValid()) { + pvs.emplace_back(pv); + } + } + } + return pvs; + } +} diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ.cc deleted file mode 100644 index 956a3eab1eb48..0000000000000 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ.cc +++ /dev/null @@ -1,713 +0,0 @@ -#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h" -#include "RecoVertex/VertexPrimitives/interface/VertexException.h" - -using namespace std; - -namespace { - - bool recTrackLessZ1(const DAClusterizerInZ::track_t& tk1, const DAClusterizerInZ::track_t& tk2) { - return tk1.z < tk2.z; - } -} // namespace - -vector DAClusterizerInZ::fill(const vector& tracks) const { - // prepare track data for clustering - vector tks; - for (vector::const_iterator it = tracks.begin(); it != tracks.end(); it++) { - track_t t; - t.z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z(); - double tantheta = tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta()); - double phi = ((*it).stateAtBeamLine().trackStateAtPCA()).momentum().phi(); - // get the beam-spot - reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot(); - t.dz2 = pow((*it).track().dzError(), 2) // track errror - + (pow(beamspot.BeamWidthX() * cos(phi), 2) + pow(beamspot.BeamWidthY() * sin(phi), 2)) / - pow(tantheta, 2) // beam-width induced - + pow(vertexSize_, 2); // intrinsic vertex size, safer for outliers and short lived decays - if (d0CutOff_ > 0) { - Measurement1D IP = (*it).stateAtBeamLine().transverseImpactParameter(); // error constains beamspot - t.pi = 1. / (1. + exp(pow(IP.value() / IP.error(), 2) - pow(d0CutOff_, 2))); // reduce weight for high ip tracks - } else { - t.pi = 1.; - } - t.tt = &(*it); - t.Z = 1.; - tks.push_back(t); - } - return tks; -} - -double DAClusterizerInZ::Eik(const track_t& t, const vertex_t& k) const { return pow(t.z - k.z, 2) / t.dz2; } - -double DAClusterizerInZ::update(double beta, vector& tks, vector& y) const { - //update weights and vertex positions - // mass constrained annealing without noise - // returns the squared sum of changes of vertex positions - - unsigned int nt = tks.size(); - - //initialize sums - double sumpi = 0; - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->se = 0; - k->sw = 0; - k->swz = 0; - k->swE = 0; - k->Tc = 0; - } - - // loop over tracks - for (unsigned int i = 0; i < nt; i++) { - // update pik and Zi - double Zi = 0.; - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->ei = exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time - Zi += k->pk * k->ei; - } - tks[i].Z = Zi; - - // normalization for pk - if (tks[i].Z > 0) { - sumpi += tks[i].pi; - // accumulate weighted z and weights for vertex update - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->se += tks[i].pi * k->ei / Zi; - double w = k->pk * tks[i].pi * k->ei / Zi / tks[i].dz2; - k->sw += w; - k->swz += w * tks[i].z; - k->swE += w * Eik(tks[i], *k); - } - } else { - sumpi += tks[i].pi; - } - - } // end of track loop - - // now update z and pk - double delta = 0; - for (vector::iterator k = y.begin(); k != y.end(); k++) { - if (k->sw > 0) { - double znew = k->swz / k->sw; - delta += pow(k->z - znew, 2); - k->z = znew; - k->Tc = 2 * k->swE / k->sw; - } else { - edm::LogInfo("sumw") << "invalid sum of weights in fit: " << k->sw << endl; - if (verbose_) { - cout << " a cluster melted away ? pk=" << k->pk << " sumw=" << k->sw << endl; - } - k->Tc = -1; - } - - k->pk = k->pk * k->se / sumpi; - } - - // return how much the prototypes moved - return delta; -} - -double DAClusterizerInZ::update(double beta, vector& tks, vector& y, double& rho0) const { - // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise - // returns the squared sum of changes of vertex positions - - unsigned int nt = tks.size(); - - //initialize sums - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->se = 0; - k->sw = 0; - k->swz = 0; - k->swE = 0; - k->Tc = 0; - } - - // loop over tracks - for (unsigned int i = 0; i < nt; i++) { - // update pik and Zi - double Zi = rho0 * exp(-beta * dzCutOff_ * dzCutOff_); // cut-off - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->ei = exp(-beta * Eik(tks[i], *k)); // cache exponential for one track at a time - Zi += k->pk * k->ei; - } - tks[i].Z = Zi; - - // normalization - if (tks[i].Z > 0) { - // accumulate weighted z and weights for vertex update - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->se += tks[i].pi * k->ei / Zi; - double w = k->pk * tks[i].pi * k->ei / Zi / tks[i].dz2; - k->sw += w; - k->swz += w * tks[i].z; - k->swE += w * Eik(tks[i], *k); - } - } - - } // end of track loop - - // now update z - double delta = 0; - for (vector::iterator k = y.begin(); k != y.end(); k++) { - if (k->sw > 0) { - double znew = k->swz / k->sw; - delta += pow(k->z - znew, 2); - k->z = znew; - k->Tc = 2 * k->swE / k->sw; - } else { - edm::LogInfo("sumw") << "invalid sum of weights in fit: " << k->sw << endl; - if (verbose_) { - cout << " a cluster melted away ? pk=" << k->pk << " sumw=" << k->sw << endl; - } - k->Tc = 0; - } - } - - // return how much the prototypes moved - return delta; -} - -bool DAClusterizerInZ::merge(vector& y, int nt) const { - // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise - - if (y.size() < 2) - return false; - - for (vector::iterator k = y.begin(); (k + 1) != y.end(); k++) { - if (fabs((k + 1)->z - k->z) < 1.e-3) { // with fabs if only called after freeze-out (splitAll() at highter T) - double rho = k->pk + (k + 1)->pk; - if (rho > 0) { - k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho; - } else { - k->z = 0.5 * (k->z + (k + 1)->z); - } - k->pk = rho; - - y.erase(k + 1); - return true; - } - } - - return false; -} - -bool DAClusterizerInZ::merge(vector& y, double& beta) const { - // merge clusters that collapsed or never separated, - // only merge if the estimated critical temperature of the merged vertex is below the current temperature - // return true if vertices were merged, false otherwise - if (y.size() < 2) - return false; - - for (vector::iterator k = y.begin(); (k + 1) != y.end(); k++) { - if (fabs((k + 1)->z - k->z) < 2.e-3) { - double rho = k->pk + (k + 1)->pk; - double swE = k->swE + (k + 1)->swE - k->pk * (k + 1)->pk / rho * pow((k + 1)->z - k->z, 2); - double Tc = 2 * swE / (k->sw + (k + 1)->sw); - - if (Tc * beta < 1) { - if (rho > 0) { - k->z = (k->pk * k->z + (k + 1)->z * (k + 1)->pk) / rho; - } else { - k->z = 0.5 * (k->z + (k + 1)->z); - } - k->pk = rho; - k->sw += (k + 1)->sw; - k->swE = swE; - k->Tc = Tc; - y.erase(k + 1); - return true; - } - } - } - - return false; -} - -bool DAClusterizerInZ::purge(vector& y, vector& tks, double& rho0, const double beta) const { - // eliminate clusters with only one significant/unique track - if (y.size() < 2) - return false; - - unsigned int nt = tks.size(); - double sumpmin = nt; - vector::iterator k0 = y.end(); - for (vector::iterator k = y.begin(); k != y.end(); k++) { - int nUnique = 0; - double sump = 0; - double pmax = k->pk / (k->pk + rho0 * exp(-beta * dzCutOff_ * dzCutOff_)); - for (unsigned int i = 0; i < nt; i++) { - if (tks[i].Z > 0) { - double p = k->pk * exp(-beta * Eik(tks[i], *k)) / tks[i].Z; - sump += p; - if ((p > 0.9 * pmax) && (tks[i].pi > 0)) { - nUnique++; - } - } - } - - if ((nUnique < 2) && (sump < sumpmin)) { - sumpmin = sump; - k0 = k; - } - } - - if (k0 != y.end()) { - if (verbose_) { - cout << "eliminating prototype at " << k0->z << " with sump=" << sumpmin << endl; - } - //rho0+=k0->pk; - y.erase(k0); - return true; - } else { - return false; - } -} - -double DAClusterizerInZ::beta0(double betamax, vector& tks, vector& y) const { - double T0 = 0; // max Tc for beta=0 - // estimate critical temperature from beta=0 (T=inf) - unsigned int nt = tks.size(); - - for (vector::iterator k = y.begin(); k != y.end(); k++) { - // vertex fit at T=inf - double sumwz = 0; - double sumw = 0; - for (unsigned int i = 0; i < nt; i++) { - double w = tks[i].pi / tks[i].dz2; - sumwz += w * tks[i].z; - sumw += w; - } - k->z = sumwz / sumw; - - // estimate Tcrit, eventually do this in the same loop - double a = 0, b = 0; - for (unsigned int i = 0; i < nt; i++) { - double dx = tks[i].z - (k->z); - double w = tks[i].pi / tks[i].dz2; - a += w * pow(dx, 2) / tks[i].dz2; - b += w; - } - double Tc = 2. * a / b; // the critical temperature of this vertex - if (Tc > T0) - T0 = Tc; - } // vertex loop (normally there should be only one vertex at beta=0) - - if (T0 > 1. / betamax) { - return betamax / pow(coolingFactor_, int(log(T0 * betamax) / log(coolingFactor_)) - 1); - } else { - // ensure at least one annealing step - return betamax / coolingFactor_; - } -} - -bool DAClusterizerInZ::split(double beta, vector& tks, vector& y, double threshold) const { - // split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1) - // an update must have been made just before doing this (same beta, no merging) - // returns true if at least one cluster was split - - double epsilon = 1e-3; // split all single vertices by 10 um - bool split = false; - - // avoid left-right biases by splitting highest Tc first - - std::vector > critical; - for (unsigned int ik = 0; ik < y.size(); ik++) { - if (beta * y[ik].Tc > 1.) { - critical.push_back(make_pair(y[ik].Tc, ik)); - } - } - stable_sort(critical.begin(), critical.end(), std::greater >()); - - for (unsigned int ic = 0; ic < critical.size(); ic++) { - unsigned int ik = critical[ic].second; - // estimate subcluster positions and weight - double p1 = 0, z1 = 0, w1 = 0; - double p2 = 0, z2 = 0, w2 = 0; - //double sumpi=0; - for (unsigned int i = 0; i < tks.size(); i++) { - if (tks[i].Z > 0) { - //sumpi+=tks[i].pi; - double p = y[ik].pk * exp(-beta * Eik(tks[i], y[ik])) / tks[i].Z * tks[i].pi; - double w = p / tks[i].dz2; - if (tks[i].z < y[ik].z) { - p1 += p; - z1 += w * tks[i].z; - w1 += w; - } else { - p2 += p; - z2 += w * tks[i].z; - w2 += w; - } - } - } - if (w1 > 0) { - z1 = z1 / w1; - } else { - z1 = y[ik].z - epsilon; - } - if (w2 > 0) { - z2 = z2 / w2; - } else { - z2 = y[ik].z + epsilon; - } - - // reduce split size if there is not enough room - if ((ik > 0) && (y[ik - 1].z >= z1)) { - z1 = 0.5 * (y[ik].z + y[ik - 1].z); - } - if ((ik + 1 < y.size()) && (y[ik + 1].z <= z2)) { - z2 = 0.5 * (y[ik].z + y[ik + 1].z); - } - - // split if the new subclusters are significantly separated - if ((z2 - z1) > epsilon) { - split = true; - vertex_t vnew; - vnew.pk = p1 * y[ik].pk / (p1 + p2); - y[ik].pk = p2 * y[ik].pk / (p1 + p2); - vnew.z = z1; - y[ik].z = z2; - y.insert(y.begin() + ik, vnew); - - // adjust remaining pointers - for (unsigned int jc = ic; jc < critical.size(); jc++) { - if (critical[jc].second > ik) { - critical[jc].second++; - } - } - } - } - - // stable_sort(y.begin(), y.end(), clusterLessZ); - return split; -} - -void DAClusterizerInZ::splitAll(vector& y) const { - double epsilon = 1e-3; // split all single vertices by 10 um - double zsep = 2 * epsilon; // split vertices that are isolated by at least zsep (vertices that haven't collapsed) - vector y1; - - for (vector::iterator k = y.begin(); k != y.end(); k++) { - if (((k == y.begin()) || (k - 1)->z < k->z - zsep) && (((k + 1) == y.end()) || (k + 1)->z > k->z + zsep)) { - // isolated prototype, split - vertex_t vnew; - vnew.z = k->z - epsilon; - (*k).z = k->z + epsilon; - vnew.pk = 0.5 * (*k).pk; - (*k).pk = 0.5 * (*k).pk; - y1.push_back(vnew); - y1.push_back(*k); - - } else if (y1.empty() || (y1.back().z < k->z - zsep)) { - y1.push_back(*k); - } else { - y1.back().z -= epsilon; - k->z += epsilon; - y1.push_back(*k); - } - } // vertex loop - - y = y1; -} - -DAClusterizerInZ::DAClusterizerInZ(const edm::ParameterSet& conf) { - // some defaults to avoid uninitialized variables - verbose_ = conf.getUntrackedParameter("verbose", false); - useTc_ = true; - betamax_ = 0.1; - betastop_ = 1.0; - coolingFactor_ = 0.8; - maxIterations_ = 100; - vertexSize_ = 0.05; // 0.5 mm - dzCutOff_ = 4.0; // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes - - // configure - - double Tmin = conf.getParameter("Tmin"); - vertexSize_ = conf.getParameter("vertexSize"); - coolingFactor_ = conf.getParameter("coolingFactor"); - d0CutOff_ = conf.getParameter("d0CutOff"); - dzCutOff_ = conf.getParameter("dzCutOff"); - maxIterations_ = 100; - if (Tmin == 0) { - cout << "DAClusterizerInZ: invalid Tmin" << Tmin << " reset do default " << 1. / betamax_ << endl; - } else { - betamax_ = 1. / Tmin; - } - - // for testing, negative cooling factor: revert to old splitting scheme - if (coolingFactor_ < 0) { - coolingFactor_ = -coolingFactor_; - useTc_ = false; - } -} - -void DAClusterizerInZ::dump(const double beta, - const vector& y, - const vector& tks0, - int verbosity) const { - // copy and sort for nicer printout - vector tks; - for (vector::const_iterator t = tks0.begin(); t != tks0.end(); t++) { - tks.push_back(*t); - } - stable_sort(tks.begin(), tks.end(), recTrackLessZ1); - - cout << "-----DAClusterizerInZ::dump ----" << endl; - cout << "beta=" << beta << " betamax= " << betamax_ << endl; - cout << " z= "; - cout.precision(4); - for (vector::const_iterator k = y.begin(); k != y.end(); k++) { - cout << setw(8) << fixed << k->z; - } - cout << endl << "T=" << setw(15) << 1. / beta << " Tc= "; - for (vector::const_iterator k = y.begin(); k != y.end(); k++) { - cout << setw(8) << fixed << k->Tc; - } - - cout << endl << " pk="; - for (vector::const_iterator k = y.begin(); k != y.end(); k++) { - cout << setw(8) << setprecision(3) << fixed << k->pk; - } - cout << endl; - - if (verbosity > 0) { - double E = 0, F = 0; - cout << endl; - cout << "---- z +/- dz ip +/-dip pt phi eta weights ----" << endl; - cout.precision(4); - for (unsigned int i = 0; i < tks.size(); i++) { - if (tks[i].Z > 0) { - F -= log(tks[i].Z) / beta; - } - double tz = tks[i].z; - cout << setw(3) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6) << sqrt(tks[i].dz2); - - if (tks[i].tt->track().quality(reco::TrackBase::highPurity)) { - cout << " *"; - } else { - cout << " "; - } - if (tks[i].tt->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) { - cout << "+"; - } else { - cout << "-"; - } - cout << setw(1) - << tks[i] - .tt->track() - .hitPattern() - .pixelBarrelLayersWithMeasurement(); // see DataFormats/TrackReco/interface/HitPattern.h - cout << setw(1) << tks[i].tt->track().hitPattern().pixelEndcapLayersWithMeasurement(); - cout << setw(1) << hex - << tks[i].tt->track().hitPattern().trackerLayersWithMeasurement() - - tks[i].tt->track().hitPattern().pixelLayersWithMeasurement() - << dec; - cout << "=" << setw(1) << hex - << tks[i].tt->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec; - - Measurement1D IP = tks[i].tt->stateAtBeamLine().transverseImpactParameter(); - cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error(); - cout << " " << setw(6) << setprecision(2) << tks[i].tt->track().pt() * tks[i].tt->track().charge(); - cout << " " << setw(5) << setprecision(2) << tks[i].tt->track().phi() << " " << setw(5) << setprecision(2) - << tks[i].tt->track().eta(); - - for (vector::const_iterator k = y.begin(); k != y.end(); k++) { - if ((tks[i].pi > 0) && (tks[i].Z > 0)) { - //double p=pik(beta,tks[i],*k); - double p = k->pk * exp(-beta * Eik(tks[i], *k)) / tks[i].Z; - if (p > 0.0001) { - cout << setw(8) << setprecision(3) << p; - } else { - cout << " . "; - } - E += p * Eik(tks[i], *k); - } else { - cout << " "; - } - } - cout << endl; - } - cout << endl << "T=" << 1 / beta << " E=" << E << " n=" << y.size() << " F= " << F << endl << "----------" << endl; - } -} - -vector DAClusterizerInZ::vertices(const vector& tracks, - const int verbosity) const { - vector tks = fill(tracks); - unsigned int nt = tracks.size(); - double rho0 = 0.0; // start with no outlier rejection - - vector clusters; - if (tks.empty()) - return clusters; - - vector y; // the vertex prototypes - - // initialize:single vertex at infinite temperature - vertex_t vstart; - vstart.z = 0.; - vstart.pk = 1.; - y.push_back(vstart); - int niter = 0; // number of iterations - - // estimate first critical temperature - double beta = beta0(betamax_, tks, y); - niter = 0; - while ((update(beta, tks, y) > 1.e-6) && (niter++ < maxIterations_)) { - } - - // annealing loop, stop when T1/Tmin) - while (beta < betamax_) { - if (useTc_) { - update(beta, tks, y); - while (merge(y, beta)) { - update(beta, tks, y); - } - split(beta, tks, y, 1.); - beta = beta / coolingFactor_; - } else { - beta = beta / coolingFactor_; - splitAll(y); - } - - // make sure we are not too far from equilibrium before cooling further - niter = 0; - while ((update(beta, tks, y) > 1.e-6) && (niter++ < maxIterations_)) { - } - } - - if (useTc_) { - // last round of splitting, make sure no critical clusters are left - update(beta, tks, y); - while (merge(y, beta)) { - update(beta, tks, y); - } - unsigned int ntry = 0; - while (split(beta, tks, y, 1.) && (ntry++ < 10)) { - niter = 0; - while ((update(beta, tks, y) > 1.e-6) && (niter++ < maxIterations_)) { - } - merge(y, beta); - update(beta, tks, y); - } - } else { - // merge collapsed clusters - while (merge(y, beta)) { - update(beta, tks, y); - } - if (verbose_) { - cout << "dump after 1st merging " << endl; - dump(beta, y, tks, 2); - } - } - - // switch on outlier rejection - rho0 = 1. / nt; - for (vector::iterator k = y.begin(); k != y.end(); k++) { - k->pk = 1.; - } // democratic - niter = 0; - while ((update(beta, tks, y, rho0) > 1.e-8) && (niter++ < maxIterations_)) { - } - if (verbose_) { - cout << "rho0=" << rho0 << " niter=" << niter << endl; - dump(beta, y, tks, 2); - } - - // merge again (some cluster split by outliers collapse here) - while (merge(y, tks.size())) { - } - if (verbose_) { - cout << "dump after 2nd merging " << endl; - dump(beta, y, tks, 2); - } - - // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices - while (beta <= betastop_) { - while (purge(y, tks, rho0, beta)) { - niter = 0; - while ((update(beta, tks, y, rho0) > 1.e-6) && (niter++ < maxIterations_)) { - } - } - beta /= coolingFactor_; - niter = 0; - while ((update(beta, tks, y, rho0) > 1.e-6) && (niter++ < maxIterations_)) { - } - } - - // // new, one last round of cleaning at T=Tstop - // while(purge(y,tks,rho0, beta)){ - // niter=0; while((update(beta, tks,y,rho0) > 1.e-6) && (niter++ < maxIterations_)){ } - // } - - if (verbose_) { - cout << "Final result, rho0=" << rho0 << endl; - dump(beta, y, tks, 2); - } - - // select significant tracks and use a TransientVertex as a container - GlobalError dummyError; - - // ensure correct normalization of probabilities, should make double assginment reasonably impossible - for (unsigned int i = 0; i < nt; i++) { - tks[i].Z = rho0 * exp(-beta * dzCutOff_ * dzCutOff_); - for (vector::iterator k = y.begin(); k != y.end(); k++) { - tks[i].Z += k->pk * exp(-beta * Eik(tks[i], *k)); - } - } - - for (vector::iterator k = y.begin(); k != y.end(); k++) { - GlobalPoint pos(0, 0, k->z); - vector vertexTracks; - for (unsigned int i = 0; i < nt; i++) { - if (tks[i].Z > 0) { - double p = k->pk * exp(-beta * Eik(tks[i], *k)) / tks[i].Z; - if ((tks[i].pi > 0) && (p > 0.5)) { - vertexTracks.push_back(*(tks[i].tt)); - tks[i].Z = 0; - } // setting Z=0 excludes double assignment - } - } - TransientVertex v(pos, dummyError, vertexTracks, 5); - clusters.push_back(v); - } - - return clusters; -} - -vector > DAClusterizerInZ::clusterize(const vector& tracks) const { - if (verbose_) { - cout << "###################################################" << endl; - cout << "# DAClusterizerInZ::clusterize nt=" << tracks.size() << endl; - cout << "###################################################" << endl; - } - - vector > clusters; - vector pv = vertices(tracks); - - if (verbose_) { - cout << "# DAClusterizerInZ::clusterize pv.size=" << pv.size() << endl; - } - if (pv.empty()) { - return clusters; - } - - // fill into clusters and merge - vector aCluster = pv.begin()->originalTracks(); - - for (vector::iterator k = pv.begin() + 1; k != pv.end(); k++) { - if (fabs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) { - // close a cluster - clusters.push_back(aCluster); - aCluster.clear(); - } - for (unsigned int i = 0; i < k->originalTracks().size(); i++) { - aCluster.push_back(k->originalTracks().at(i)); - } - } - clusters.push_back(aCluster); - - return clusters; -} diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc index 7e96ab91f9c09..61ef41b672e80 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZT_vect.cc @@ -238,7 +238,7 @@ DAClusterizerInZT_vect::track_t DAClusterizerInZT_vect::fill(const vector 0.3) || (std::abs(t_t) > t0Max_)) { + if ((tk.dtErrorExt() > TransientTrackBuilder::defaultInvalidTrackTimeReso) || (std::abs(t_t) > t0Max_)) { t_dt2 = 0; // tracks with no time measurement } else { t_dt2 = 1. / t_dt2; @@ -1560,7 +1560,7 @@ void DAClusterizerInZT_vect::dump(const double beta, const vertex_t& y, const tr void DAClusterizerInZT_vect::fillPSetDescription(edm::ParameterSetDescription& desc) { DAClusterizerInZ_vect::fillPSetDescription(desc); - desc.add("tmerge", 0.01); // 4D only + desc.add("tmerge", 0.1); // 4D only desc.add("dtCutOff", 4.); // 4D only desc.add("t0Max", 1.0); // 4D only desc.add("vertexSizeTime", 0.008); // 4D only diff --git a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc index 72470403eee5e..acefd0b8dd200 100644 --- a/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc +++ b/RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc @@ -48,7 +48,7 @@ DAClusterizerInZ_vect::DAClusterizerInZ_vect(const edm::ParameterSet& conf) { overlap_frac_ = conf.getParameter("overlap_frac"); #ifdef DEBUG - std::cout << "DAClusterizerinZ_vect: mintrkweight = " << mintrkweight << std::endl; + std::cout << "DAClusterizerinZ_vect: mintrkweight = " << mintrkweight_ << std::endl; std::cout << "DAClusterizerinZ_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl; std::cout << "DAClusterizerInZ_vect: uniquetrkminp = " << uniquetrkminp_ << std::endl; std::cout << "DAClusterizerinZ_vect: zmerge = " << zmerge_ << std::endl; @@ -63,6 +63,10 @@ DAClusterizerInZ_vect::DAClusterizerInZ_vect(const edm::ParameterSet& conf) { std::cout << "DAClusterizerinZ_vect: convergence mode = " << convergence_mode_ << std::endl; std::cout << "DAClusterizerinZ_vect: delta_highT = " << delta_highT_ << std::endl; std::cout << "DAClusterizerinZ_vect: delta_lowT = " << delta_lowT_ << std::endl; + + std::cout << "DAClusterizerinZ_vect: run in blocks = " << runInBlocks_ << std::endl; + std::cout << "DAClusterizerinZ_vect: block_size = " << block_size_ << std::endl; + std::cout << "DAClusterizerinZ_vect: overlap_fraction = " << overlap_frac_ << std::endl; std::cout << "DAClusterizerinZ_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl; #endif @@ -212,7 +216,6 @@ DAClusterizerInZ_vect::track_t DAClusterizerInZ_vect::fill(const vector 1) { std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.zvtx[k0] - << " with sump=" << sumpmin << " rho*nt =" << y.rho[k0] * nt << endl; + << " with sump=" << sumpmin << " rho*nt =" << y.rho[k0] * nt << " pnUnique=" << pnUnique[k0] << endl; } #endif @@ -763,11 +766,10 @@ bool DAClusterizerInZ_vect::split(const double beta, track_t& tks, vertex_t& y, return split; } -vector DAClusterizerInZ_vect::vertices(const vector& tracks) const { +vector DAClusterizerInZ_vect::vertices_no_blocks(const vector& tracks) const { track_t&& tks = fill(tracks); tks.extractRaw(); - unsigned int nt = tks.getSize(); double rho0 = 0.0; // start with no outlier rejection vector clusters; @@ -807,7 +809,7 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "last round of splitting" << std::endl; } #endif @@ -835,7 +837,7 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "turning on outlier rejection at T=" << 1 / beta << std::endl; } #endif @@ -853,7 +855,7 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "merging with outlier rejection at T=" << 1 / beta << std::endl; } if (DEBUGLEVEL > 2) @@ -869,7 +871,7 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "after merging with outlier rejection at T=" << 1 / beta << std::endl; } if (DEBUGLEVEL > 2) @@ -898,7 +900,7 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "last cooling T=" << 1 / beta << std::endl; } #endif @@ -912,66 +914,15 @@ vector DAClusterizerInZ_vect::vertices(const vector 0) { - std::cout << "DAClusterizerInZ_vect::vertices :" + std::cout << "DAClusterizerInZ_vect::vertices_no_blocks :" << "stop cooling at T=" << 1 / beta << std::endl; } if (DEBUGLEVEL > 2) dump(beta, y, tks, 2, rho0); #endif - // select significant tracks and use a TransientVertex as a container - - set_vtx_range(beta, tks, y); - const unsigned int nv = y.getSize(); - for (unsigned int k = 0; k < nv; k++) { - if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) { - y.rho[k] = 0; - y.zvtx[k] = 0; - } - } - - const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); - std::vector> vtx_track_indices(nv); - for (unsigned int i = 0; i < nt; i++) { - const auto kmin = tks.kmin[i]; - const auto kmax = tks.kmax[i]; - for (auto k = kmin; k < kmax; k++) { - y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i]); - } - - local_exp_list_range(y.exp_arg, y.exp, kmin, kmax); - - tks.sum_Z[i] = z_sum_init; - for (auto k = kmin; k < kmax; k++) { - tks.sum_Z[i] += y.rho[k] * y.exp[k]; - } - const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0; - - for (auto k = kmin; k < kmax; k++) { - double p = y.rho[k] * y.exp[k] * invZ; - if (p > mintrkweight_) { - // assign track i -> vertex k (hard, mintrkweight should be >= 0.5 here - vtx_track_indices[k].push_back(i); - break; - } - } - - } // track loop - - GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); - for (unsigned int k = 0; k < nv; k++) { - if (!vtx_track_indices[k].empty()) { - GlobalPoint pos(0, 0, y.zvtx[k]); - vector vertexTracks; - for (auto i : vtx_track_indices[k]) { - vertexTracks.push_back(*(tks.tt[i])); - } - TransientVertex v(pos, dummyError, vertexTracks, 0); - clusters.push_back(v); - } - } - - return clusters; + // assign tracks and fill into transient vertices + return fill_vertices(beta, rho0, tks, y); } vector DAClusterizerInZ_vect::vertices_in_blocks(const vector& tracks) const { @@ -1241,8 +1192,9 @@ vector DAClusterizerInZ_vect::vertices_in_blocks(const vector DAClusterizerInZ_vect::vertices_in_blocks(const vector clusters; + if (nv == 0) { + return clusters; + } GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); + vector vertexTracks; + for (unsigned int k = 0; k < nv; k++) { if (!vtx_track_indices[k].empty()) { - GlobalPoint pos(0, 0, vertices_tot[k].first); - vector vertexTracks; for (auto i : vtx_track_indices[k]) { vertexTracks.push_back(*(tracks_tot.tt[i])); #ifdef DEBUG - std::cout << y.zvtx[k] << "," << (*tks.tt[i]).stateAtBeamLine().trackStateAtPCA().position().z() << std::endl; + std::cout << vertices_tot[k].first << "," + << (*(tracks_tot.tt[i])).stateAtBeamLine().trackStateAtPCA().position().z() << std::endl; #endif } - TransientVertex v(pos, dummyError, vertexTracks, 0); - clusters.push_back(v); + } + + // implement what clusterize() did before : merge left-to-right if distance < 2 * vertexSize_ + if ((k + 1 == nv) || (abs(vertices_tot[k + 1].first - vertices_tot[k].first) > (2 * vertexSize_))) { + // close a cluster + if (vertexTracks.size() > 1) { + GlobalPoint pos(0, 0, vertices_tot[k].first); // only usable with subsequent fit + TransientVertex v(pos, dummyError, vertexTracks, 0); + clusters.push_back(v); + } + vertexTracks.clear(); } } return clusters; -} +} // end of vertices_in_blocks -vector> DAClusterizerInZ_vect::clusterize( - const vector& tracks) const { - vector> clusters; +vector DAClusterizerInZ_vect::fill_vertices(double beta, double rho0, track_t& tks, vertex_t& y) const { + // select significant tracks and use a TransientVertex as a container - if (tracks.empty()) - return clusters; + set_vtx_range(beta, tks, y); + const unsigned int nv = y.getSize(); + for (unsigned int k = 0; k < nv; k++) { + if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) { + y.rho[k] = 0; + y.zvtx[k] = 0; + } + } - vector pv; + // ensure consistent assignment probabillities and make a hard assignment + const unsigned int nt = tks.getSize(); + const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_); + std::vector> vtx_track_indices(nv); + std::vector> vtx_track_weights(nv); + for (unsigned int i = 0; i < nt; i++) { + const auto kmin = tks.kmin[i]; + const auto kmax = tks.kmax[i]; + for (auto k = kmin; k < kmax; k++) { + y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i]); + } + + local_exp_list_range(y.exp_arg, y.exp, kmin, kmax); + + tks.sum_Z[i] = z_sum_init; + for (auto k = kmin; k < kmax; k++) { + tks.sum_Z[i] += y.rho[k] * y.exp[k]; + } + const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0; + + double pmax = -1; + unsigned int k_pmax = 0; + for (auto k = kmin; k < kmax; k++) { + double p = y.rho[k] * y.exp[k] * invZ; + if (p > pmax) { + pmax = p; + k_pmax = k; + } + } + + if (pmax > mintrkweight_) { + // assign to the cluster with the highest assignment weight, if it is at least mintrkweight_ + vtx_track_indices[k_pmax].push_back(i); + vtx_track_weights[k_pmax].push_back(pmax); + } + } + + // fill transient vertices + // the position is normally not used, probably not optimal when Tstop <> 2, anyway + vector clusters; + for (unsigned int k = 0; k < nv; k++) { + double sump = 0; + double sumw = 0; + double sumwp = 0, sumwz = 0; + if (!vtx_track_indices[k].empty()) { + vector vertexTracks; + TransientVertex::TransientTrackToFloatMap trkWeightMap; + unsigned int j = 0; + for (auto i : vtx_track_indices[k]) { + auto p = vtx_track_weights[k][j]; + vertexTracks.push_back(*(tks.tt[i])); + trkWeightMap[vertexTracks[j]] = p; + auto w = p * tks.dz2[i]; + sump += p; + sumw += w; + sumwp += w * p; + sumwz += w * tks.zpca[i]; + j++; + } + float zerror_squared = 1.; // + if ((sumw > 0) && (sumwp > 0)) { + zerror_squared = sumwp / (sumw * sumw); + y.zvtx[k] = sumwz / sumw; + } + const auto& bs = vertexTracks[0].stateAtBeamLine().beamSpot(); + GlobalPoint pos(bs.x(y.zvtx[k]), bs.y(y.zvtx[k]), y.zvtx[k]); + const float xerror_squared = pow(bs.BeamWidthX(), 2); + const float yerror_squared = pow(bs.BeamWidthY(), 2); + GlobalError err(xerror_squared, 0, yerror_squared, 0., 0., zerror_squared); + TransientVertex v(pos, err, vertexTracks, 0, 2 * sump - 3.); + v.weightMap(trkWeightMap); + clusters.push_back(v); + } + } + + return clusters; +} + +vector DAClusterizerInZ_vect::vertices(const vector& tracks) const { if (runInBlocks_ and (block_size_ < tracks.size())) //doesn't bother if low number of tracks - pv = vertices_in_blocks(tracks); + return vertices_in_blocks(tracks); else - pv = vertices(tracks); + return vertices_no_blocks(tracks); +} + +vector> DAClusterizerInZ_vect::clusterize( // OBSOLETE + const vector& tracks) const { + vector> clusters; + vector&& pv = vertices(tracks); #ifdef DEBUG if (DEBUGLEVEL > 0) { @@ -1494,6 +1548,6 @@ void DAClusterizerInZ_vect::fillPSetDescription(edm::ParameterSetDescription& de desc.add("uniquetrkminp", 0.0); desc.add("zrange", 4.0); desc.add("runInBlocks", false); - desc.add("block_size", 512); - desc.add("overlap_frac", 0.5); + desc.add("block_size", 10000); + desc.add("overlap_frac", 0.0); } diff --git a/RecoVertex/PrimaryVertexProducer/src/GapClusterizerInZ.cc b/RecoVertex/PrimaryVertexProducer/src/GapClusterizerInZ.cc index 9a3d6c0b33ed9..495718ea25dee 100644 --- a/RecoVertex/PrimaryVertexProducer/src/GapClusterizerInZ.cc +++ b/RecoVertex/PrimaryVertexProducer/src/GapClusterizerInZ.cc @@ -61,6 +61,20 @@ vector > GapClusterizerInZ::clusterize(const vector return clusters; } +vector GapClusterizerInZ::vertices(const vector& tracks) const { + /* repackage track clusters, compatibility with newer clusterizers */ + std::vector primary_vertices; + auto trackClusters = clusterize(tracks); + + GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01); + for (auto& vertexTracks : trackClusters) { + GlobalPoint position(0, 0, 0); // dummy + primary_vertices.push_back(TransientVertex(position, dummyError, vertexTracks, 0)); + } + + return primary_vertices; +} + void GapClusterizerInZ::fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add("zSeparation", 1.0); desc.addUntracked("verbose", false); diff --git a/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc b/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc index 70d502aae3de4..74ae7c36154ce 100644 --- a/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc +++ b/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc @@ -39,9 +39,6 @@ PrimaryVertexProducerAlgorithm::PrimaryVertexProducerAlgorithm(const edm::Parame if (clusteringAlgorithm == "gap") { theTrackClusterizer = new GapClusterizerInZ( conf.getParameter("TkClusParameters").getParameter("TkGapClusParameters")); - } else if (clusteringAlgorithm == "DA") { - theTrackClusterizer = new DAClusterizerInZ( - conf.getParameter("TkClusParameters").getParameter("TkDAClusParameters")); } // provide the vectorized version of the clusterizer, if supported by the build else if (clusteringAlgorithm == "DA_vect") { diff --git a/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmFromTracksPID.cc b/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmFromTracksPID.cc new file mode 100644 index 0000000000000..db1fd074b4621 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmFromTracksPID.cc @@ -0,0 +1,193 @@ +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ValidatedPluginMacros.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "vdt/vdtMath.h" + +#include "RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h" + +#ifdef PVTX_DEBUG +#define LOG edm::LogPrint("VertexTimeAlgorithmFromTracksPID") +#else +#define LOG LogDebug("VertexTimeAlgorithmFromTracksPID") +#endif + +VertexTimeAlgorithmFromTracksPID::VertexTimeAlgorithmFromTracksPID(edm::ParameterSet const& iConfig, + edm::ConsumesCollector& iCC) + : VertexTimeAlgorithmBase(iConfig, iCC), + trackMTDTimeToken_(iCC.consumes(iConfig.getParameter("trackMTDTimeVMapTag"))), + trackMTDTimeErrorToken_(iCC.consumes(iConfig.getParameter("trackMTDTimeErrorVMapTag"))), + trackMTDTimeQualityToken_(iCC.consumes(iConfig.getParameter("trackMTDTimeQualityVMapTag"))), + trackMTDTofPiToken_(iCC.consumes(iConfig.getParameter("trackMTDTofPiVMapTag"))), + trackMTDTofKToken_(iCC.consumes(iConfig.getParameter("trackMTDTofKVMapTag"))), + trackMTDTofPToken_(iCC.consumes(iConfig.getParameter("trackMTDTofPVMapTag"))), + minTrackVtxWeight_(iConfig.getParameter("minTrackVtxWeight")), + minTrackTimeQuality_(iConfig.getParameter("minTrackTimeQuality")), + probPion_(iConfig.getParameter("probPion")), + probKaon_(iConfig.getParameter("probKaon")), + probProton_(iConfig.getParameter("probProton")), + Tstart_(iConfig.getParameter("Tstart")), + coolingFactor_(iConfig.getParameter("coolingFactor")) {} + +void VertexTimeAlgorithmFromTracksPID::fillPSetDescription(edm::ParameterSetDescription& iDesc) { + VertexTimeAlgorithmBase::fillPSetDescription(iDesc); + + iDesc.add("trackMTDTimeVMapTag", edm::InputTag("trackExtenderWithMTD:generalTracktmtd")) + ->setComment("Input ValueMap for track time at MTD"); + iDesc.add("trackMTDTimeErrorVMapTag", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd")) + ->setComment("Input ValueMap for track time uncertainty at MTD"); + iDesc.add("trackMTDTimeQualityVMapTag", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA")) + ->setComment("Input ValueMap for track MVA quality value"); + iDesc.add("trackMTDTofPiVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi")) + ->setComment("Input ValueMap for track tof as pion"); + iDesc.add("trackMTDTofKVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackTofK")) + ->setComment("Input ValueMap for track tof as kaon"); + iDesc.add("trackMTDTofPVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackTofP")) + ->setComment("Input ValueMap for track tof as proton"); + + iDesc.add("minTrackVtxWeight", 0.5)->setComment("Minimum track weight"); + iDesc.add("minTrackTimeQuality", 0.8)->setComment("Minimum MVA Quality selection on tracks"); + + iDesc.add("probPion", 0.7)->setComment("A priori probability pions"); + iDesc.add("probKaon", 0.2)->setComment("A priori probability kaons"); + iDesc.add("probProton", 0.1)->setComment("A priori probability protons"); + + iDesc.add("Tstart", 256.)->setComment("DA initial temperature T"); + iDesc.add("coolingFactor", 0.5)->setComment("DA cooling factor"); +} + +void VertexTimeAlgorithmFromTracksPID::setEvent(edm::Event& iEvent, edm::EventSetup const&) { + // additional collections required for vertex-time calculation + trackMTDTimes_ = iEvent.get(trackMTDTimeToken_); + trackMTDTimeErrors_ = iEvent.get(trackMTDTimeErrorToken_); + trackMTDTimeQualities_ = iEvent.get(trackMTDTimeQualityToken_); + trackMTDTofPi_ = iEvent.get(trackMTDTofPiToken_); + trackMTDTofK_ = iEvent.get(trackMTDTofKToken_); + trackMTDTofP_ = iEvent.get(trackMTDTofPToken_); +} + +bool VertexTimeAlgorithmFromTracksPID::vertexTime(float& vtxTime, + float& vtxTimeError, + const TransientVertex& vtx) const { + if (vtx.originalTracks().empty()) { + return false; + } + + auto const vtxTime_init = vtxTime; + auto const vtxTimeError_init = vtxTimeError; + const int max_iterations = 100; + + double tsum = 0; + double wsum = 0; + double w2sum = 0; + + double const a[3] = {probPion_, probKaon_, probProton_}; + + std::vector v_trackInfo; + v_trackInfo.reserve(vtx.originalTracks().size()); + + // initial guess + for (const auto& trk : vtx.originalTracks()) { + auto const trkWeight = vtx.trackWeight(trk); + if (trkWeight > minTrackVtxWeight_) { + auto const trkTimeQuality = trackMTDTimeQualities_[trk.trackBaseRef()]; + + if (trkTimeQuality >= minTrackTimeQuality_) { + auto const trkTime = trackMTDTimes_[trk.trackBaseRef()]; + auto const trkTimeError = trackMTDTimeErrors_[trk.trackBaseRef()]; + + v_trackInfo.emplace_back(); + auto& trkInfo = v_trackInfo.back(); + + trkInfo.trkWeight = trkWeight; + trkInfo.trkTimeError = trkTimeError; + + trkInfo.trkTimeHyp[0] = trkTime - trackMTDTofPi_[trk.trackBaseRef()]; + trkInfo.trkTimeHyp[1] = trkTime - trackMTDTofK_[trk.trackBaseRef()]; + trkInfo.trkTimeHyp[2] = trkTime - trackMTDTofP_[trk.trackBaseRef()]; + + auto const wgt = trkWeight / (trkTimeError * trkTimeError); + wsum += wgt; + + for (uint j = 0; j < 3; ++j) { + tsum += wgt * trkInfo.trkTimeHyp[j] * a[j]; + } + LOG << "vertexTimeFromTracks: track" + << " pt=" << trk.track().pt() << " eta=" << trk.track().eta() << " phi=" << trk.track().phi() + << " vtxWeight=" << trkWeight << " time=" << trkTime << " timeError=" << trkTimeError + << " timeQuality=" << trkTimeQuality << " timeHyp[pion]=" << trkInfo.trkTimeHyp[0] + << " timeHyp[kaon]=" << trkInfo.trkTimeHyp[1] << " timeHyp[proton]=" << trkInfo.trkTimeHyp[2]; + } + } + } + if (wsum > 0) { + auto t0 = tsum / wsum; + auto beta = 1. / Tstart_; + int nit = 0; + while ((nit++) < max_iterations) { + tsum = 0; + wsum = 0; + w2sum = 0; + + for (auto const& trkInfo : v_trackInfo) { + double dt = trkInfo.trkTimeError; + double e[3] = {0, 0, 0}; + const double cut_off = 4.5; + double Z = vdt::fast_exp( + -beta * cut_off); // outlier rejection term Z_0 = exp(-beta * cut_off) = exp(-beta * 0.5 * 3 * 3) + for (unsigned int j = 0; j < 3; j++) { + auto const tpull = (trkInfo.trkTimeHyp[j] - t0) / dt; + e[j] = vdt::fast_exp(-0.5 * beta * tpull * tpull); + Z += a[j] * e[j]; + } + + double wsum_trk = 0; + for (uint j = 0; j < 3; j++) { + double wt = a[j] * e[j] / Z; + double w = wt * trkInfo.trkWeight / (dt * dt); + wsum_trk += w; + tsum += w * trkInfo.trkTimeHyp[j]; + } + + wsum += wsum_trk; + w2sum += wsum_trk * wsum_trk * (dt * dt) / trkInfo.trkWeight; + } + + if (wsum < 1e-10) { + LOG << "vertexTimeFromTracks: failed while iterating"; + return false; + } + + vtxTime = tsum / wsum; + + LOG << "vertexTimeFromTracks: iteration=" << nit << ", T= " << 1 / beta << ", t=" << vtxTime + << ", t-t0=" << vtxTime - t0; + + if ((std::abs(vtxTime - t0) < 1e-4 / std::sqrt(beta)) and beta >= 1.) { + vtxTimeError = std::sqrt(w2sum) / wsum; + + LOG << "vertexTimeFromTracks: tfit = " << vtxTime << " +/- " << vtxTimeError << " trec = " << vtx.time() + << ", iteration=" << nit; + + return true; + } + + if ((std::abs(vtxTime - t0) < 1e-3) and beta < 1.) { + beta = std::min(1., beta / coolingFactor_); + } + + t0 = vtxTime; + } + + LOG << "vertexTimeFromTracks: failed to converge"; + } else { + LOG << "vertexTimeFromTracks: has no track timing info"; + } + + vtxTime = vtxTime_init; + vtxTimeError = vtxTimeError_init; + + return false; +} diff --git a/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmLegacy4D.cc b/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmLegacy4D.cc new file mode 100644 index 0000000000000..b6dd746fecf77 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmLegacy4D.cc @@ -0,0 +1,58 @@ +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmLegacy4D.h" + +#ifdef PVTX_DEBUG +#define LOG edm::LogPrint("VertexTimeAlgorithmLegacy4D") +#else +#define LOG LogDebug("VertexTimeAlgorithmLegacy4D") +#endif + +VertexTimeAlgorithmLegacy4D::VertexTimeAlgorithmLegacy4D(edm::ParameterSet const& iConfig, edm::ConsumesCollector& iCC) + : VertexTimeAlgorithmBase(iConfig, iCC) {} + +void VertexTimeAlgorithmLegacy4D::fillPSetDescription(edm::ParameterSetDescription& iDesc) { + VertexTimeAlgorithmBase::fillPSetDescription(iDesc); +} + +void VertexTimeAlgorithmLegacy4D::setEvent(edm::Event& iEvent, edm::EventSetup const&){}; + +bool VertexTimeAlgorithmLegacy4D::vertexTime(float& vtxTime, float& vtxTimeError, const TransientVertex& vtx) const { + const auto num_track = vtx.originalTracks().size(); + if (num_track == 0) { + return false; + } + + double sumwt = 0.; + double sumwt2 = 0.; + double sumw = 0.; + double vartime = 0.; + + for (const auto& trk : vtx.originalTracks()) { + const double time = trk.timeExt(); + const double err = trk.dtErrorExt(); + if ((time == 0) && (err > TransientTrackBuilder::defaultInvalidTrackTimeReso)) + continue; // tracks with no time information, as implemented in TransientTrackBuilder.cc l.17 + const double inverr = err > 0. ? 1.0 / err : 0.; + const double w = inverr * inverr; + sumwt += w * time; + sumwt2 += w * time * time; + sumw += w; + } + + if (sumw > 0) { + double sumsq = sumwt2 - sumwt * sumwt / sumw; + double chisq = num_track > 1 ? sumsq / double(num_track - 1) : sumsq / double(num_track); + vartime = chisq / sumw; + + vtxTime = sumwt / sumw; + vtxTimeError = sqrt(vartime); + return true; + } + + vtxTime = 0; + vtxTimeError = 1.; + return false; +} diff --git a/TrackingTools/TransientTrack/interface/TransientTrackBuilder.h b/TrackingTools/TransientTrack/interface/TransientTrackBuilder.h index 21280d9c2d5e4..8e20fffcc995d 100644 --- a/TrackingTools/TransientTrack/interface/TransientTrackBuilder.h +++ b/TrackingTools/TransientTrack/interface/TransientTrackBuilder.h @@ -69,6 +69,7 @@ class TransientTrackBuilder { const MagneticField* field() const { return theField; } const edm::ESHandle trackingGeometry() const { return theTrackingGeometry; } + static constexpr float defaultInvalidTrackTimeReso = 0.350f; private: const MagneticField* theField; diff --git a/TrackingTools/TransientTrack/src/TransientTrackBuilder.cc b/TrackingTools/TransientTrack/src/TransientTrackBuilder.cc index afe55fa17f3dc..748b15814938c 100644 --- a/TrackingTools/TransientTrack/src/TransientTrackBuilder.cc +++ b/TrackingTools/TransientTrack/src/TransientTrackBuilder.cc @@ -13,10 +13,6 @@ using namespace reco; using namespace std; using namespace edm; -namespace { - constexpr float defaultInvalidTrackReso = 0.350f; -} - TransientTrack TransientTrackBuilder::build(const Track* t) const { return TransientTrack(*t, theField, theTrackingGeometry); } @@ -105,11 +101,11 @@ vector TransientTrackBuilder::build(const edm::Handle 1e-6 ? timeReso : defaultInvalidTrackReso); // make the error much larger than the BS time width + timeReso = (timeReso > 1e-6 ? timeReso + : defaultInvalidTrackTimeReso); // make the error much larger than the BS time width if (edm::isNotFinite(time)) { time = 0.0; - timeReso = defaultInvalidTrackReso; + timeReso = defaultInvalidTrackTimeReso; } ttVect.push_back(TransientTrack(ref, time, timeReso, theField, theTrackingGeometry)); } @@ -125,11 +121,11 @@ vector TransientTrackBuilder::build(const edm::Handle 1e-6 ? timeReso : defaultInvalidTrackReso); // make the error much larger than the BS time width + timeReso = (timeReso > 1e-6 ? timeReso + : defaultInvalidTrackTimeReso); // make the error much larger than the BS time width if (edm::isNotFinite(time)) { time = 0.0; - timeReso = defaultInvalidTrackReso; + timeReso = defaultInvalidTrackTimeReso; } ttVect.push_back(TransientTrack(new GsfTransientTrack(ref, time, timeReso, theField, theTrackingGeometry))); } @@ -148,11 +144,11 @@ vector TransientTrackBuilder::build(const edm::Handle(trkColl, i).castTo(); double time = trackTimes[ref]; double timeReso = trackTimeResos[ref]; - timeReso = - (timeReso > 1e-6 ? timeReso : defaultInvalidTrackReso); // make the error much larger than the BS time width + timeReso = (timeReso > 1e-6 ? timeReso + : defaultInvalidTrackTimeReso); // make the error much larger than the BS time width if (edm::isNotFinite(time)) { time = 0.0; - timeReso = defaultInvalidTrackReso; + timeReso = defaultInvalidTrackTimeReso; } ttVect.push_back(TransientTrack(new GsfTransientTrack( RefToBase(trkColl, i).castTo(), time, timeReso, theField, theTrackingGeometry))); @@ -160,11 +156,11 @@ vector TransientTrackBuilder::build(const edm::Handle(trkColl, i).castTo(); double time = trackTimes[ref]; double timeReso = trackTimeResos[ref]; - timeReso = - (timeReso > 1e-6 ? timeReso : defaultInvalidTrackReso); // make the error much larger than the BS time width + timeReso = (timeReso > 1e-6 ? timeReso + : defaultInvalidTrackTimeReso); // make the error much larger than the BS time width if (edm::isNotFinite(time)) { time = 0.0; - timeReso = defaultInvalidTrackReso; + timeReso = defaultInvalidTrackTimeReso; } ttVect.push_back(TransientTrack( RefToBase(trkColl, i).castTo(), time, timeReso, theField, theTrackingGeometry)); diff --git a/Validation/Configuration/python/valid3Dt.py b/Validation/Configuration/python/valid3Dt.py new file mode 100644 index 0000000000000..0d15cb8771d83 --- /dev/null +++ b/Validation/Configuration/python/valid3Dt.py @@ -0,0 +1,17 @@ +import FWCore.ParameterSet.Config as cms + +def customise(process): + process.mtdTracksValid.inputTagV = 'offlinePrimaryVertices' + process.mtdTracksValid.t0SafePID = 'tofPID3D:t0safe' + process.mtdTracksValid.sigmat0SafePID = 'tofPID3D:sigmat0safe' + process.mtdTracksValid.sigmat0PID = 'tofPID3D:sigmat0' + process.mtdTracksValid.t0PID = 'tofPID3D:t0' + + process.vertices4DValid.offline4DPV = 'offlinePrimaryVertices' + process.vertices4DValid.t0PID = 'tofPID3D:t0' + process.vertices4DValid.t0SafePID = 'tofPID3D:t0safe' + process.vertices4DValid.sigmat0SafePID = 'tofPID3D:sigmat0safe' + process.vertices4DValid.probPi = 'tofPID3D:probPi' + process.vertices4DValid.probK = 'tofPID3D:probK' + process.vertices4DValid.probP = 'tofPID3D:probP' + return(process)