From 6278d0078b2ff350bae24971c9efeed6aa9a8a87 Mon Sep 17 00:00:00 2001 From: abaty Date: Sat, 24 Nov 2018 19:59:16 +0100 Subject: [PATCH] Adding Vertex Recovery for peripheral events for 2018 PbPb (#174) * Vertex Recovery added for 2018 PbPb data * updated other configs * moved checking for vtx before transientTrack builder to increase speed significantly * changed from vString to vInputTag for mass replace to work --- .../test/runForestAOD_pponAA_DATA_103X.py | 9 + .../test/runForestAOD_pponAA_MB_103X.py | 8 + .../test/runForestAOD_pponAA_MIX_103X.py | 7 + .../TrackAnalysis/python/TrkAnalyzers_cff.py | 2 +- .../TrackAnalysis/python/trackAnalyzer_cfi.py | 2 +- .../TrackAnalysis/src/TrackAnalyzer.cc | 2 +- .../interface/PrimaryVertexRecoveryProducer.h | 95 +++++ .../plugins/PrimaryVertexRecoveryProducer.cc | 380 ++++++++++++++++++ .../OfflinePrimaryVerticesRecovery_cfi.py | 49 +++ 9 files changed, 551 insertions(+), 3 deletions(-) create mode 100644 RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexRecoveryProducer.h create mode 100644 RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexRecoveryProducer.cc create mode 100644 RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVerticesRecovery_cfi.py diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_DATA_103X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_DATA_103X.py index 081e70ade9628..d9bc949f531a0 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_DATA_103X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_DATA_103X.py @@ -166,12 +166,17 @@ process.rechitanalyzerpp.zdcRecHitSrc = cms.untracked.InputTag("QWzdcreco") ############################################################################### +#Recover peripheral primary vertices +#https://twiki.cern.ch/twiki/bin/view/CMS/HITracking2018PbPb#Peripheral%20Vertex%20Recovery +process.load("RecoVertex.PrimaryVertexProducer.OfflinePrimaryVerticesRecovery_cfi") + ######################### # Main analysis list ######################### process.ana_step = cms.Path( + process.offlinePrimaryVerticesRecovery + process.HiForest + process.hltanalysis + process.hltobject + @@ -248,6 +253,10 @@ process.pAna = cms.EndPath(process.skimanalysis) +from HLTrigger.Configuration.CustomConfigs import MassReplaceInputTag +process = MassReplaceInputTag(process,"offlinePrimaryVertices","offlinePrimaryVerticesRecovery") +process.offlinePrimaryVerticesRecovery.oldVertexLabel = "offlinePrimaryVertices" + ############################################################################### # Customization diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MB_103X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MB_103X.py index d46105867436c..1c23aeb13af48 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MB_103X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MB_103X.py @@ -169,12 +169,16 @@ process.load('HeavyIonsAnalysis.JetAnalysis.rechitanalyzer_cfi') ############################################################################### +#Recover peripheral primary vertices +#https://twiki.cern.ch/twiki/bin/view/CMS/HITracking2018PbPb#Peripheral%20Vertex%20Recovery +process.load("RecoVertex.PrimaryVertexProducer.OfflinePrimaryVerticesRecovery_cfi") ######################### # Main analysis list ######################### process.ana_step = cms.Path( + process.offlinePrimaryVerticesRecovery + process.HiForest + process.runAnalyzer + process.hltanalysis + @@ -251,5 +255,9 @@ process.pAna = cms.EndPath(process.skimanalysis) +from HLTrigger.Configuration.CustomConfigs import MassReplaceInputTag +process = MassReplaceInputTag(process,"offlinePrimaryVertices","offlinePrimaryVerticesRecovery") +process.offlinePrimaryVerticesRecovery.oldVertexLabel = "offlinePrimaryVertices" + # Customization ############################################################################### diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MIX_103X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MIX_103X.py index d0af564685efb..2a1831a5fdf7d 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MIX_103X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pponAA_MIX_103X.py @@ -167,12 +167,16 @@ process.load('HeavyIonsAnalysis.JetAnalysis.rechitanalyzer_cfi') ############################################################################### +#Recover peripheral primary vertices +#https://twiki.cern.ch/twiki/bin/view/CMS/HITracking2018PbPb#Peripheral%20Vertex%20Recovery +process.load("RecoVertex.PrimaryVertexProducer.OfflinePrimaryVerticesRecovery_cfi") ######################### # Main analysis list ######################### process.ana_step = cms.Path( + process.offlinePrimaryVerticesRecovery + process.HiForest + process.runAnalyzer + process.hltanalysis + @@ -249,5 +253,8 @@ process.pAna = cms.EndPath(process.skimanalysis) +from HLTrigger.Configuration.CustomConfigs import MassReplaceInputTag +process = MassReplaceInputTag(process,"offlinePrimaryVertices","offlinePrimaryVerticesRecovery") +process.offlinePrimaryVerticesRecovery.oldVertexLabel = "offlinePrimaryVertices" # Customization ############################################################################### diff --git a/HeavyIonsAnalysis/TrackAnalysis/python/TrkAnalyzers_cff.py b/HeavyIonsAnalysis/TrackAnalysis/python/TrkAnalyzers_cff.py index 42a4bdb6fee60..cdc9fa0431083 100644 --- a/HeavyIonsAnalysis/TrackAnalysis/python/TrkAnalyzers_cff.py +++ b/HeavyIonsAnalysis/TrackAnalysis/python/TrkAnalyzers_cff.py @@ -5,7 +5,7 @@ anaTrack = ppTrack.clone( trackPtMin = 0.49, trackSrc = cms.InputTag("hiGeneralTracks"), - vertexSrc = cms.vstring('hiSelectedVertex'), + vertexSrc = cms.VInputTag('hiSelectedVertex'), mvaSrc = cms.InputTag('hiGeneralTracks','MVAVals'), pfCandSrc = cms.InputTag("particleFlowTmp"), doMVA = False diff --git a/HeavyIonsAnalysis/TrackAnalysis/python/trackAnalyzer_cfi.py b/HeavyIonsAnalysis/TrackAnalysis/python/trackAnalyzer_cfi.py index e7233a9086f1a..a0d780e3635a5 100644 --- a/HeavyIonsAnalysis/TrackAnalysis/python/trackAnalyzer_cfi.py +++ b/HeavyIonsAnalysis/TrackAnalysis/python/trackAnalyzer_cfi.py @@ -4,7 +4,7 @@ 'TrackAnalyzer', trackPtMin = cms.untracked.double(0.01), simTrackPtMin = cms.untracked.double(0.1), - vertexSrc = cms.vstring('offlinePrimaryVertices'), + vertexSrc = cms.VInputTag('offlinePrimaryVertices'), trackSrc = cms.InputTag('generalTracks'), mvaSrc = cms.InputTag("generalTracks","MVAValues"), particleSrc = cms.InputTag('genParticles'), diff --git a/HeavyIonsAnalysis/TrackAnalysis/src/TrackAnalyzer.cc b/HeavyIonsAnalysis/TrackAnalysis/src/TrackAnalyzer.cc index 0cda2589f7017..6fcb85216b396 100644 --- a/HeavyIonsAnalysis/TrackAnalysis/src/TrackAnalyzer.cc +++ b/HeavyIonsAnalysis/TrackAnalysis/src/TrackAnalyzer.cc @@ -343,7 +343,7 @@ TrackAnalyzer::TrackAnalyzer(const edm::ParameterSet& iConfig) associatorMapRS_ = consumes(iConfig.getParameter("associatorMap")); } - std::vector vertexSrcString_ = iConfig.getParameter>("vertexSrc"); + std::vector vertexSrcString_ = iConfig.getParameter>("vertexSrc"); for (const auto& src : vertexSrcString_) { vertexSrc_.push_back(consumes(edm::InputTag(src))); } diff --git a/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexRecoveryProducer.h b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexRecoveryProducer.h new file mode 100644 index 0000000000000..68319d393cf8f --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexRecoveryProducer.h @@ -0,0 +1,95 @@ +// -*- C++ -*- +// +// Package: PrimaryVertexProducer +// Class: PrimaryVertexProducer +// +/**\class PrimaryVertexProducer PrimaryVertexProducer.cc RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducer.cc + + Description: steers tracker primary vertex reconstruction and storage + + Implementation: + +*/ +// +// Original Author: Pascal Vanlaer +// Created: Tue Feb 28 11:06:34 CET 2006 +// +// + + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +//#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFindingBase.h" +#include "RecoVertex/PrimaryVertexProducer/interface/TrackClusterizerInZ.h" +#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h" +#include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZT_vect.h" + + +#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" +#include "RecoVertex/VertexTools/interface/VertexCompatibleWithBeam.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +// +// class declaration +// + +class PrimaryVertexRecoveryProducer : public edm::stream::EDProducer<> { +public: + PrimaryVertexRecoveryProducer(const edm::ParameterSet&); + ~PrimaryVertexRecoveryProducer() override; + + void produce(edm::Event&, const edm::EventSetup&) override; + + // access to config + edm::ParameterSet config() const { return theConfig; } + +private: + // ----------member data --------------------------- + TrackFilterForPVFindingBase* theTrackFilter; + TrackClusterizerInZ* theTrackClusterizer; + + // vtx fitting algorithms + struct algo { + VertexFitter<5> * fitter; + VertexCompatibleWithBeam * vertexSelector; + std::string label; + bool useBeamConstraint; + double minNdof; + }; + + std::vector< algo > algorithms; + + edm::ParameterSet theConfig; + bool fVerbose; + + edm::EDGetTokenT bsToken; + bool redoAllVertices; + edm::EDGetTokenT oldVtxToken; + edm::EDGetTokenT trkToken; + edm::EDGetTokenT > trkTimesToken; + edm::EDGetTokenT > trkTimeResosToken; + + bool f4D; +}; diff --git a/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexRecoveryProducer.cc b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexRecoveryProducer.cc new file mode 100644 index 0000000000000..d6c40ca9d7827 --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/plugins/PrimaryVertexRecoveryProducer.cc @@ -0,0 +1,380 @@ +#include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexRecoveryProducer.h" + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" + +#include "RecoVertex/VertexTools/interface/GeometricAnnealing.h" +#include + +PrimaryVertexRecoveryProducer::PrimaryVertexRecoveryProducer(const edm::ParameterSet& conf) + :theConfig(conf) +{ + + fVerbose = conf.getUntrackedParameter("verbose", false); + + redoAllVertices = conf.getParameter("redoAllVertices"); + oldVtxToken = consumes(conf.getParameter("oldVertexLabel")); + + trkToken = consumes(conf.getParameter("TrackLabel")); + bsToken = consumes(conf.getParameter("beamSpotLabel")); + f4D = false; + + // select and configure the track selection + std::string trackSelectionAlgorithm=conf.getParameter("TkFilterParameters").getParameter("algorithm"); + if(trackSelectionAlgorithm=="filter"){ + theTrackFilter= new TrackFilterForPVFinding( conf.getParameter("TkFilterParameters") ); + }else if (trackSelectionAlgorithm=="filterWithThreshold"){ + theTrackFilter= new HITrackFilterForPVFinding(conf.getParameter("TkFilterParameters")); + }else{ + throw VertexException("PrimaryVertexRecoveryProducerAlgorithm: unknown track selection algorithm: " + trackSelectionAlgorithm); + } + + + // select and configure the track clusterizer + std::string clusteringAlgorithm=conf.getParameter("TkClusParameters").getParameter("algorithm"); + 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") { + 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{ + throw VertexException("PrimaryVertexRecoveryProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm); + } + + if( f4D ) { + trkTimesToken = consumes >(conf.getParameter("TrackTimesLabel") ); + trkTimeResosToken = consumes >(conf.getParameter("TrackTimeResosLabel") ); + } + + + // select and configure the vertex fitters + if (conf.exists("vertexCollections")){ + std::vector vertexCollections =conf.getParameter< std::vector >("vertexCollections"); + + for( std::vector< edm::ParameterSet >::const_iterator algoconf = vertexCollections.begin(); algoconf != vertexCollections.end(); algoconf++){ + + algo algorithm; + std::string fitterAlgorithm = algoconf->getParameter("algorithm"); + if (fitterAlgorithm=="KalmanVertexFitter") { + algorithm.fitter= new KalmanVertexFitter(); + } else if( fitterAlgorithm=="AdaptiveVertexFitter") { + algorithm.fitter= new AdaptiveVertexFitter( GeometricAnnealing( algoconf->getParameter("chi2cutoff"))); + } else { + throw VertexException("PrimaryVertexRecoveryProducerAlgorithm: 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")); + algorithms.push_back(algorithm); + + produces(algorithm.label); + } + }else{ + edm::LogWarning("MisConfiguration")<<"this module's configuration has changed, please update to have a vertexCollections=cms.VPSet parameter."; + + algo algorithm; + std::string fitterAlgorithm = conf.getParameter("algorithm"); + if (fitterAlgorithm=="KalmanVertexFitter") { + algorithm.fitter= new KalmanVertexFitter(); + } else if( fitterAlgorithm=="AdaptiveVertexFitter") { + algorithm.fitter= new AdaptiveVertexFitter(); + } else { + throw VertexException("PrimaryVertexRecoveryProducerAlgorithm: unknown algorithm: " + fitterAlgorithm); + } + algorithm.label = ""; + algorithm.minNdof = conf.getParameter("minNdof"); + algorithm.useBeamConstraint=conf.getParameter("useBeamConstraint"); + + algorithm.vertexSelector=new VertexCompatibleWithBeam(VertexDistanceXY(), conf.getParameter("PVSelParameters").getParameter("maxDistanceToBeam")); + + algorithms.push_back(algorithm); + produces(algorithm.label); + } + + +} + + +PrimaryVertexRecoveryProducer::~PrimaryVertexRecoveryProducer() { + if (theTrackFilter) delete theTrackFilter; + if (theTrackClusterizer) delete theTrackClusterizer; + for( std::vector ::const_iterator algorithm=algorithms.begin(); algorithm!=algorithms.end(); algorithm++){ + if (algorithm->fitter) delete algorithm->fitter; + if (algorithm->vertexSelector) delete algorithm->vertexSelector; + } +} + + +void +PrimaryVertexRecoveryProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + + // get the BeamSpot, it will alwys be needed, even when not used as a constraint + reco::BeamSpot beamSpot; + edm::Handle recoBeamSpotHandle; + iEvent.getByToken(bsToken,recoBeamSpotHandle); + if (recoBeamSpotHandle.isValid()){ + beamSpot = *recoBeamSpotHandle; + }else{ + edm::LogError("UnusableBeamSpot") << "No beam spot available from EventSetup"; + } + + bool validBS = true; + 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 "< tks; + iEvent.getByToken(trkToken, tks); + + //check if we alreayd have a vertex + unsigned int nAlgosFound = 0; + for( std::vector ::const_iterator algorithm=algorithms.begin(); algorithm!=algorithms.end(); algorithm++){ + auto result = std::make_unique(); + reco::VertexCollection & vColl = (*result); + + if(!redoAllVertices){ + //see if this event already has a vertex + const reco::VertexCollection* recoVertices; + edm::Handle< reco::VertexCollection > vtxs; + iEvent.getByToken(oldVtxToken, vtxs); + recoVertices = vtxs.product(); + if(recoVertices->size()!=0){ + bool hasValid = false; + unsigned int maxNtrk = 0; + unsigned int biggestIndx = 0; + for(size_t i = 0; i < recoVertices->size(); ++i){ + if(!((*recoVertices)[i].isFake())){ + hasValid = true; + if((*recoVertices)[i].nTracks() > maxNtrk){ + maxNtrk = (*recoVertices)[i].nTracks(); + biggestIndx = i; + } + } + } + if(hasValid){ + nAlgosFound++; + vColl.push_back((*recoVertices)[biggestIndx]); + iEvent.put(std::move(result), algorithm->label); + continue; + } + } + } + } + if(nAlgosFound==algorithms.size()) return; + + // interface RECO tracks to vertex reconstruction + edm::ESHandle theB; + iSetup.get().get("TransientTrackBuilder",theB); + 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())); + } else { + t_tks = (*theB).build(tks, beamSpot); + } + if(fVerbose) {std::cout << "RecoVertex/PrimaryVertexRecoveryProducer" + << "Found: " << t_tks.size() << " reconstructed tracks" << "\n"; + } + + + // select tracks + std::vector && seltks = theTrackFilter->select( t_tks ); + + // clusterize tracks in Z + std::vector< std::vector > && clusters = theTrackClusterizer->clusterize(seltks); + + if (fVerbose){std::cout << " clustering returned "<< clusters.size() << " clusters from " << seltks.size() << " selected tracks" <::const_iterator algorithm=algorithms.begin(); algorithm!=algorithms.end(); algorithm++){ + + + auto result = std::make_unique(); + reco::VertexCollection & vColl = (*result); + + std::vector pvs; + for (std::vector< std::vector >::const_iterator iclus + = clusters.begin(); iclus != clusters.end(); iclus++) { + + double meantime = 0.; + double expv_x2 = 0.; + double normw = 0.; + if( f4D ) { + for( const auto& tk : *iclus ) { + const double time = tk.timeExt(); + const double inverr = 1.0/tk.dtErrorExt(); + const double w = inverr*inverr; + meantime += time*w; + expv_x2 += time*time*w; + normw += w; + } + meantime = meantime/normw; + expv_x2 = expv_x2/normw; + } + const double time_var = ( f4D ? expv_x2 - meantime*meantime : 0. ); + + + TransientVertex v; + if( algorithm->useBeamConstraint && validBS &&((*iclus).size()>1) ){ + + v = algorithm->fitter->vertex(*iclus, beamSpot); + + if( f4D ) { + if( v.isValid() ) { + auto err = v.positionError().matrix4D(); + err(3,3) = time_var/(double)iclus->size(); + v = TransientVertex(v.position(),meantime,err,v.originalTracks(),v.totalChiSquared()); + } + } + + }else if( !(algorithm->useBeamConstraint) && ((*iclus).size()>1) ) { + + v = algorithm->fitter->vertex(*iclus); + + if( f4D ) { + if( v.isValid() ) { + auto err = v.positionError().matrix4D(); + err(3,3) = time_var/(double)iclus->size(); + v = TransientVertex(v.position(),meantime,err,v.originalTracks(),v.totalChiSquared()); + } + } + + }// else: no fit ==> v.isValid()=False + + + 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; + } + } + + if ( v.isValid() + && (v.degreesOfFreedom()>=algorithm->minNdof) + && (!validBS || (*(algorithm->vertexSelector))(v,beamVertexState)) + ) pvs.push_back(v); + }// end of cluster loop + + if(fVerbose){ + std::cout << "PrimaryVertexRecoveryProducerAlgorithm::vertices candidates =" << pvs.size() << std::endl; + } + + + if (clusters.size()>2 && clusters.size() > 2*pvs.size()) + edm::LogWarning("PrimaryVertexRecoveryProducer") << "more than half of candidate vertices lost " << pvs.size() << ' ' << clusters.size(); + + if (pvs.empty() && seltks.size()>5) + edm::LogWarning("PrimaryVertexRecoveryProducer") << "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()); + } + + + + // 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 (vColl.empty()) { + GlobalError bse(beamSpot.rotatedCovariance3D()); + if ( (bse.cxx() <= 0.) || + (bse.cyy() <= 0.) || + (bse.czz() <= 0.) ) { + AlgebraicSymMatrix33 we; + we(0,0)=10000; we(1,1)=10000; we(2,2)=10000; + vColl.push_back(reco::Vertex(beamSpot.position(), we,0.,0.,0)); + if(fVerbose){ + std::cout <<"RecoVertex/PrimaryVertexRecoveryProducer: " + << "Beamspot with invalid errors "<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(); + } + std::cout << std::endl; + } + } + + if(!(vColl[0].isFake())){ + std::cout << "Recovered a Vertex with " << std::setw(3) << vColl[0].tracksSize() << " tracks!" << std::endl; + std::cout << "x: " << std::setw(6) << vColl[0].position().x() << " y: " << std::setw(6) << vColl[0].position().y() << " z: " << std::setw(6) << vColl[0].position().z() << " chi2/ndof: " << std::setw(4) << vColl[0].chi2()/vColl[0].ndof() << std::endl; + } + iEvent.put(std::move(result), algorithm->label); + } + +} + + +//define this as a plug-in +DEFINE_FWK_MODULE(PrimaryVertexRecoveryProducer); diff --git a/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVerticesRecovery_cfi.py b/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVerticesRecovery_cfi.py new file mode 100644 index 0000000000000..31973d0ec951f --- /dev/null +++ b/RecoVertex/PrimaryVertexProducer/python/OfflinePrimaryVerticesRecovery_cfi.py @@ -0,0 +1,49 @@ +import FWCore.ParameterSet.Config as cms + +offlinePrimaryVerticesRecovery = cms.EDProducer( + "PrimaryVertexRecoveryProducer", + + redoAllVertices = cms.bool(False),#Only turn this on if you want to spend lots of CPU time, nt recommended + oldVertexLabel = cms.InputTag("offlinePrimaryVertices"), + + verbose = cms.untracked.bool(False), + TrackLabel = cms.InputTag("generalTracks"), + beamSpotLabel = cms.InputTag("offlineBeamSpot"), + + TkFilterParameters = cms.PSet( + algorithm=cms.string('filter'), + maxNormalizedChi2 = cms.double(999.0), + minPixelLayersWithHits=cms.int32(0), + minSiliconLayersWithHits = cms.int32(0), + maxD0Significance = cms.double(999.0), + minPt = cms.double(0.0), + maxEta = cms.double(999.0), + trackQuality = cms.string("any") + ), + + TkClusParameters = cms.PSet( + algorithm = cms.string("gap"), + TkGapClusParameters = cms.PSet( + zSeparation = cms.double(1.0) + ) + ), + + 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) + ), + #do not run vertexing with beamspot constraint + #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) + # ) + ] + ) +)