#include <RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducer.cc>

Inheritance diagram for PrimaryVertexProducer:

Classes
struct	algo

Public Member Functions
edm::ParameterSet	config () const

	PrimaryVertexProducer (const edm::ParameterSet &)

void	produce (edm::Event &, const edm::EventSetup &) override

	~PrimaryVertexProducer () override

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

	EDProducer (const EDProducer &)=delete

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () const final

const EDProducer &	operator= (const EDProducer &)=delete

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

Private Attributes
std::vector< algo >	algorithms

edm::EDGetTokenT< reco::BeamSpot >	bsToken

bool	f4D

bool	fRecoveryIteration

bool	fVerbose

edm::EDGetTokenT< reco::VertexCollection >	recoveryVtxToken

edm::ParameterSet	theConfig

TrackClusterizerInZ *	theTrackClusterizer

TrackFilterForPVFindingBase *	theTrackFilter

const edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord >	theTTBToken

edm::EDGetTokenT< edm::ValueMap< float > >	trkTimeResosToken

edm::EDGetTokenT< edm::ValueMap< float > >	trkTimesToken

edm::EDGetTokenT< reco::TrackCollection >	trkToken

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
using	CacheTypes = CacheContexts< T... >

using	GlobalCache = typename CacheTypes::GlobalCache

using	HasAbility = AbilityChecker< T... >

using	InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache

using	LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache

using	LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >

using	LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache

using	RunCache = typename CacheTypes::RunCache

using	RunContext = RunContextT< RunCache, GlobalCache >

using	RunSummaryCache = typename CacheTypes::RunSummaryCache

Detailed Description

Description: steers tracker primary vertex reconstruction and storage

Implementation: <Notes on="" implementation>="">

Definition at line 56 of file PrimaryVertexProducer.h.

Constructor & Destructor Documentation

◆ PrimaryVertexProducer()

PrimaryVertexProducer::PrimaryVertexProducer ( const edm::ParameterSet & conf )

Definition at line 18 of file PrimaryVertexProducer.cc.

     : theTTBToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), theConfig(conf) {
   fVerbose = conf.getUntrackedParameter<bool>("verbose", false);
  
   trkToken = consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("TrackLabel"));
   bsToken = consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotLabel"));
   f4D = false;
  
   // select and configure the track selection
   std::string trackSelectionAlgorithm =
       conf.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
   if (trackSelectionAlgorithm == "filter") {
     theTrackFilter = new TrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
   } else if (trackSelectionAlgorithm == "filterWithThreshold") {
     theTrackFilter = new HITrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
   } else {
     throw VertexException("PrimaryVertexProducer: unknown track selection algorithm: " + trackSelectionAlgorithm);
   }
  
   // select and configure the track clusterizer
   std::string clusteringAlgorithm =
       conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
   if (clusteringAlgorithm == "gap") {
     theTrackClusterizer = new GapClusterizerInZ(
         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
   } else if (clusteringAlgorithm == "DA") {
     theTrackClusterizer = new DAClusterizerInZ(
         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
   }
   // provide the vectorized version of the clusterizer, if supported by the build
   else if (clusteringAlgorithm == "DA_vect") {
     theTrackClusterizer = new DAClusterizerInZ_vect(
         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
   } else if (clusteringAlgorithm == "DA2D_vect") {
     theTrackClusterizer = new DAClusterizerInZT_vect(
         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
     f4D = true;
   }
  
   else {
     throw VertexException("PrimaryVertexProducer: unknown clustering algorithm: " + clusteringAlgorithm);
   }
  
   if (f4D) {
     trkTimesToken = consumes<edm::ValueMap<float> >(conf.getParameter<edm::InputTag>("TrackTimesLabel"));
     trkTimeResosToken = consumes<edm::ValueMap<float> >(conf.getParameter<edm::InputTag>("TrackTimeResosLabel"));
   }
  
   // select and configure the vertex fitters
   if (conf.exists("vertexCollections")) {
     std::vector<edm::ParameterSet> vertexCollections =
         conf.getParameter<std::vector<edm::ParameterSet> >("vertexCollections");
  
     for (std::vector<edm::ParameterSet>::const_iterator algoconf = vertexCollections.begin();
          algoconf != vertexCollections.end();
          algoconf++) {
       algo algorithm;
       std::string fitterAlgorithm = algoconf->getParameter<std::string>("algorithm");
       if (fitterAlgorithm == "KalmanVertexFitter") {
         algorithm.fitter = new KalmanVertexFitter();
       } else if (fitterAlgorithm == "AdaptiveVertexFitter") {
         algorithm.fitter = new AdaptiveVertexFitter(GeometricAnnealing(algoconf->getParameter<double>("chi2cutoff")));
       } else {
         throw VertexException("PrimaryVertexProducer: unknown algorithm: " + fitterAlgorithm);
       }
       algorithm.label = algoconf->getParameter<std::string>("label");
       algorithm.minNdof = algoconf->getParameter<double>("minNdof");
       algorithm.useBeamConstraint = algoconf->getParameter<bool>("useBeamConstraint");
       algorithm.vertexSelector =
           new VertexCompatibleWithBeam(VertexDistanceXY(), algoconf->getParameter<double>("maxDistanceToBeam"));
       algorithms.push_back(algorithm);
  
       produces<reco::VertexCollection>(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<std::string>("algorithm");
     if (fitterAlgorithm == "KalmanVertexFitter") {
       algorithm.fitter = new KalmanVertexFitter();
     } else if (fitterAlgorithm == "AdaptiveVertexFitter") {
       algorithm.fitter = new AdaptiveVertexFitter();
     } else {
       throw VertexException("PrimaryVertexProducerAlgorithm: unknown algorithm: " + fitterAlgorithm);
     }
     algorithm.label = "";
     algorithm.minNdof = conf.getParameter<double>("minNdof");
     algorithm.useBeamConstraint = conf.getParameter<bool>("useBeamConstraint");
  
     algorithm.vertexSelector = new VertexCompatibleWithBeam(
         VertexDistanceXY(),
         conf.getParameter<edm::ParameterSet>("PVSelParameters").getParameter<double>("maxDistanceToBeam"));
  
     algorithms.push_back(algorithm);
     produces<reco::VertexCollection>(algorithm.label);
   }
  
   //check if this is a recovery iteration
   fRecoveryIteration = conf.getParameter<bool>("isRecoveryIteration");
   if (fRecoveryIteration) {
     if (algorithms.empty()) {
       throw VertexException("PrimaryVertexProducer: No algorithm specified. ");
     } else if (algorithms.size() > 1) {
       throw VertexException(
           "PrimaryVertexProducer: Running in Recovery mode and more than one algorithm specified.  Please "
           "only one algorithm.");
     }
     recoveryVtxToken = consumes<reco::VertexCollection>(conf.getParameter<edm::InputTag>("recoveryVtxCollection"));
   }
 }

References qcdUeDQM_cfi::algorithm, algorithms, bsToken, edm::ParameterSet::exists(), f4D, fRecoveryIteration, fVerbose, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), recoveryVtxToken, AlCaHLTBitMon_QueryRunRegistry::string, theTrackClusterizer, theTrackFilter, trkTimeResosToken, trkTimesToken, trkToken, and HLT_FULL_cff::vertexCollections.

◆ ~PrimaryVertexProducer()

PrimaryVertexProducer::~PrimaryVertexProducer ( )

override

Definition at line 131 of file PrimaryVertexProducer.cc.

                                               {
   if (theTrackFilter)
     delete theTrackFilter;
   if (theTrackClusterizer)
     delete theTrackClusterizer;
   for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
     if (algorithm->fitter)
       delete algorithm->fitter;
     if (algorithm->vertexSelector)
       delete algorithm->vertexSelector;
   }
 }

References qcdUeDQM_cfi::algorithm, algorithms, theTrackClusterizer, and theTrackFilter.

Member Function Documentation

◆ config()

edm::ParameterSet PrimaryVertexProducer::config ( void ) const

inline

Definition at line 66 of file PrimaryVertexProducer.h.

66 { return theConfig; }

References theConfig.

◆ fillDescriptions()

void PrimaryVertexProducer::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

Definition at line 343 of file PrimaryVertexProducer.cc.

                                                                                        {
   // offlinePrimaryVertices
   edm::ParameterSetDescription desc;
   {
     edm::ParameterSetDescription vpsd1;
     vpsd1.add<double>("maxDistanceToBeam", 1.0);
     vpsd1.add<std::string>("algorithm", "AdaptiveVertexFitter");
     vpsd1.add<bool>("useBeamConstraint", false);
     vpsd1.add<std::string>("label", "");
     vpsd1.add<double>("chi2cutoff", 2.5);
     vpsd1.add<double>("minNdof", 0.0);
     std::vector<edm::ParameterSet> temp1;
     temp1.reserve(2);
     {
       edm::ParameterSet temp2;
       temp2.addParameter<double>("maxDistanceToBeam", 1.0);
       temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
       temp2.addParameter<bool>("useBeamConstraint", false);
       temp2.addParameter<std::string>("label", "");
       temp2.addParameter<double>("chi2cutoff", 2.5);
       temp2.addParameter<double>("minNdof", 0.0);
       temp1.push_back(temp2);
     }
     {
       edm::ParameterSet temp2;
       temp2.addParameter<double>("maxDistanceToBeam", 1.0);
       temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
       temp2.addParameter<bool>("useBeamConstraint", true);
       temp2.addParameter<std::string>("label", "WithBS");
       temp2.addParameter<double>("chi2cutoff", 2.5);
       temp2.addParameter<double>("minNdof", 2.0);
       temp1.push_back(temp2);
     }
     desc.addVPSet("vertexCollections", vpsd1, temp1);
   }
   desc.addUntracked<bool>("verbose", false);
   {
     edm::ParameterSetDescription psd0;
     TrackFilterForPVFinding::fillPSetDescription(psd0);
     psd0.add<int>("numTracksThreshold", 0);  // HI only
     desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
   }
   desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
   desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
   desc.add<edm::InputTag>("TrackTimeResosLabel", edm::InputTag("dummy_default"));  // 4D only
   desc.add<edm::InputTag>("TrackTimesLabel", edm::InputTag("dummy_default"));      // 4D only
  
   {
     edm::ParameterSetDescription psd0;
     {
       edm::ParameterSetDescription psd1;
       DAClusterizerInZT_vect::fillPSetDescription(psd1);
       psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
  
       edm::ParameterSetDescription psd2;
       GapClusterizerInZ::fillPSetDescription(psd2);
       psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
     }
     psd0.add<std::string>("algorithm", "DA_vect");
     desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
   }
  
   desc.add<bool>("isRecoveryIteration", false);
   desc.add<edm::InputTag>("recoveryVtxCollection", {""});
  
   descriptions.add("primaryVertexProducer", desc);
 }

References edm::ConfigurationDescriptions::add(), edm::ParameterSetDescription::add(), edm::ParameterSet::addParameter(), submitPVResolutionJobs::desc, GapClusterizerInZ::fillPSetDescription(), TrackFilterForPVFinding::fillPSetDescription(), DAClusterizerInZT_vect::fillPSetDescription(), HLT_FULL_cff::InputTag, and AlCaHLTBitMon_QueryRunRegistry::string.

◆ produce()

void PrimaryVertexProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

override

Definition at line 144 of file PrimaryVertexProducer.cc.

                                                                                  {
   // get the BeamSpot, it will always be needed, even when not used as a constraint
   reco::BeamSpot beamSpot;
   edm::Handle<reco::BeamSpot> 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 " << beamVertexState.error().matrix();
   }
  
   //if this is a recovery iteration, check if we already have a valid PV
   if (fRecoveryIteration) {
     auto const& oldVertices = iEvent.get(recoveryVtxToken);
     //look for the first valid (not-BeamSpot) vertex
     for (auto const& old : oldVertices) {
       if (!(old.isFake())) {
         //found a valid vertex, write the first one to the collection and return
         //otherwise continue with regular vertexing procedure
         auto result = std::make_unique<reco::VertexCollection>();
         result->push_back(old);
         iEvent.put(std::move(result), algorithms.begin()->label);
         return;
       }
     }
   }
  
   // get RECO tracks from the event
   // `tks` can be used as a ptr to a reco::TrackCollection
   edm::Handle<reco::TrackCollection> tks;
   iEvent.getByToken(trkToken, tks);
  
   // interface RECO tracks to vertex reconstruction
   const auto& theB = &iSetup.getData(theTTBToken);
   std::vector<reco::TransientTrack> t_tks;
  
   if (f4D) {
     edm::Handle<edm::ValueMap<float> > trackTimesH;
     edm::Handle<edm::ValueMap<float> > 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/PrimaryVertexProducer"
               << "Found: " << t_tks.size() << " reconstructed tracks"
               << "\n";
   }
  
   // select tracks
   std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
  
   // clusterize tracks in Z
   std::vector<std::vector<reco::TransientTrack> >&& clusters = theTrackClusterizer->clusterize(seltks);
  
   if (fVerbose) {
     std::cout << " clustering returned  " << clusters.size() << " clusters  from " << seltks.size()
               << " selected tracks" << std::endl;
   }
  
   // vertex fits
   for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
     auto result = std::make_unique<reco::VertexCollection>();
     reco::VertexCollection& vColl = (*result);
  
     std::vector<TransientVertex> pvs;
     for (std::vector<std::vector<reco::TransientTrack> >::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->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
  
       // 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 (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 << "PrimaryVertexProducerAlgorithm::vertices  candidates =" << pvs.size() << std::endl;
     }
  
     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());
     }
  
     // convert transient vertices returned by the theAlgo to (reco) vertices
     for (std::vector<TransientVertex>::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/PrimaryVertexProducer: "
                     << "Beamspot with invalid errors " << bse.matrix() << std::endl;
           std::cout << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n";
         }
       } 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";
         }
       }
     }
  
     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();
         }
         std::cout << std::endl;
       }
     }
  
     iEvent.put(std::move(result), algorithm->label);
   }
 }