File indexing completed on 2023-10-25 10:03:25
0001 #include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h"
0002 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0003 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008
0009 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0010 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0011 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0012
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0015
0016 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0017
0018 PrimaryVertexProducer::PrimaryVertexProducer(const edm::ParameterSet& conf)
0019 : theTTBToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), theConfig(conf) {
0020 fVerbose = conf.getUntrackedParameter<bool>("verbose", false);
0021
0022 trkToken = consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("TrackLabel"));
0023 bsToken = consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotLabel"));
0024 f4D = false;
0025 weightFit = false;
0026
0027
0028 std::string trackSelectionAlgorithm =
0029 conf.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
0030 if (trackSelectionAlgorithm == "filter") {
0031 theTrackFilter = new TrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0032 } else if (trackSelectionAlgorithm == "filterWithThreshold") {
0033 theTrackFilter = new HITrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0034 } else {
0035 throw VertexException("PrimaryVertexProducer: unknown track selection algorithm: " + trackSelectionAlgorithm);
0036 }
0037
0038
0039 std::string clusteringAlgorithm =
0040 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
0041 if (clusteringAlgorithm == "gap") {
0042 theTrackClusterizer = new GapClusterizerInZ(
0043 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
0044 } else if (clusteringAlgorithm == "DA") {
0045 theTrackClusterizer = new DAClusterizerInZ(
0046 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0047 }
0048
0049 else if (clusteringAlgorithm == "DA_vect") {
0050 theTrackClusterizer = new DAClusterizerInZ_vect(
0051 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0052 } else if (clusteringAlgorithm == "DA2D_vect") {
0053 theTrackClusterizer = new DAClusterizerInZT_vect(
0054 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0055 f4D = true;
0056 }
0057
0058 else {
0059 throw VertexException("PrimaryVertexProducer: unknown clustering algorithm: " + clusteringAlgorithm);
0060 }
0061
0062 if (f4D) {
0063 trkTimesToken = consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("TrackTimesLabel"));
0064 trkTimeResosToken = consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("TrackTimeResosLabel"));
0065 }
0066
0067
0068 std::vector<edm::ParameterSet> vertexCollections =
0069 conf.getParameter<std::vector<edm::ParameterSet>>("vertexCollections");
0070
0071 for (std::vector<edm::ParameterSet>::const_iterator algoconf = vertexCollections.begin();
0072 algoconf != vertexCollections.end();
0073 algoconf++) {
0074 algo algorithm;
0075 std::string fitterAlgorithm = algoconf->getParameter<std::string>("algorithm");
0076 if (fitterAlgorithm == "KalmanVertexFitter") {
0077 algorithm.fitter = new KalmanVertexFitter();
0078 } else if (fitterAlgorithm == "AdaptiveVertexFitter") {
0079 algorithm.fitter = new AdaptiveVertexFitter(GeometricAnnealing(algoconf->getParameter<double>("chi2cutoff")));
0080 } else if (fitterAlgorithm == "WeightedMeanFitter") {
0081 algorithm.fitter = nullptr;
0082 weightFit = true;
0083 } else {
0084 throw VertexException("PrimaryVertexProducer: unknown algorithm: " + fitterAlgorithm);
0085 }
0086 algorithm.label = algoconf->getParameter<std::string>("label");
0087 algorithm.minNdof = algoconf->getParameter<double>("minNdof");
0088 algorithm.useBeamConstraint = algoconf->getParameter<bool>("useBeamConstraint");
0089 algorithm.vertexSelector =
0090 new VertexCompatibleWithBeam(VertexDistanceXY(), algoconf->getParameter<double>("maxDistanceToBeam"));
0091 algorithms.push_back(algorithm);
0092
0093 produces<reco::VertexCollection>(algorithm.label);
0094 }
0095
0096
0097 fRecoveryIteration = conf.getParameter<bool>("isRecoveryIteration");
0098 if (fRecoveryIteration) {
0099 if (algorithms.empty()) {
0100 throw VertexException("PrimaryVertexProducer: No algorithm specified. ");
0101 } else if (algorithms.size() > 1) {
0102 throw VertexException(
0103 "PrimaryVertexProducer: Running in Recovery mode and more than one algorithm specified. Please "
0104 "only one algorithm.");
0105 }
0106 recoveryVtxToken = consumes<reco::VertexCollection>(conf.getParameter<edm::InputTag>("recoveryVtxCollection"));
0107 }
0108 }
0109
0110 PrimaryVertexProducer::~PrimaryVertexProducer() {
0111 if (theTrackFilter)
0112 delete theTrackFilter;
0113 if (theTrackClusterizer)
0114 delete theTrackClusterizer;
0115 for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0116 if (algorithm->fitter)
0117 delete algorithm->fitter;
0118 if (algorithm->vertexSelector)
0119 delete algorithm->vertexSelector;
0120 }
0121 }
0122
0123 void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0124
0125 reco::BeamSpot beamSpot;
0126 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0127 iEvent.getByToken(bsToken, recoBeamSpotHandle);
0128 if (recoBeamSpotHandle.isValid()) {
0129 beamSpot = *recoBeamSpotHandle;
0130 } else {
0131 edm::LogError("UnusableBeamSpot") << "No beam spot available from EventSetup";
0132 }
0133
0134 bool validBS = true;
0135 VertexState beamVertexState(beamSpot);
0136 if ((beamVertexState.error().cxx() <= 0.) || (beamVertexState.error().cyy() <= 0.) ||
0137 (beamVertexState.error().czz() <= 0.)) {
0138 validBS = false;
0139 edm::LogError("UnusableBeamSpot") << "Beamspot with invalid errors " << beamVertexState.error().matrix();
0140 }
0141
0142
0143 if (fRecoveryIteration) {
0144 auto const& oldVertices = iEvent.get(recoveryVtxToken);
0145
0146 for (auto const& old : oldVertices) {
0147 if (!(old.isFake())) {
0148
0149
0150 auto result = std::make_unique<reco::VertexCollection>();
0151 result->push_back(old);
0152 iEvent.put(std::move(result), algorithms.begin()->label);
0153 return;
0154 }
0155 }
0156 }
0157
0158
0159
0160 edm::Handle<reco::TrackCollection> tks;
0161 iEvent.getByToken(trkToken, tks);
0162
0163
0164 const auto& theB = &iSetup.getData(theTTBToken);
0165 std::vector<reco::TransientTrack> t_tks;
0166
0167 if (f4D) {
0168 edm::Handle<edm::ValueMap<float>> trackTimesH;
0169 edm::Handle<edm::ValueMap<float>> trackTimeResosH;
0170 iEvent.getByToken(trkTimesToken, trackTimesH);
0171 iEvent.getByToken(trkTimeResosToken, trackTimeResosH);
0172 t_tks = (*theB).build(tks, beamSpot, *(trackTimesH.product()), *(trackTimeResosH.product()));
0173 } else {
0174 t_tks = (*theB).build(tks, beamSpot);
0175 }
0176 if (fVerbose) {
0177 std::cout << "RecoVertex/PrimaryVertexProducer"
0178 << "Found: " << t_tks.size() << " reconstructed tracks"
0179 << "\n";
0180 }
0181
0182
0183 std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
0184
0185
0186 std::vector<std::vector<reco::TransientTrack>>&& clusters = theTrackClusterizer->clusterize(seltks);
0187
0188 if (fVerbose) {
0189 std::cout << " clustering returned " << clusters.size() << " clusters from " << seltks.size()
0190 << " selected tracks" << std::endl;
0191 }
0192
0193
0194 for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0195 auto result = std::make_unique<reco::VertexCollection>();
0196 reco::VertexCollection& vColl = (*result);
0197
0198 std::vector<TransientVertex> pvs;
0199 for (std::vector<std::vector<reco::TransientTrack>>::const_iterator iclus = clusters.begin();
0200 iclus != clusters.end();
0201 iclus++) {
0202 double sumwt = 0.;
0203 double sumwt2 = 0.;
0204 double sumw = 0.;
0205 double meantime = 0.;
0206 double vartime = 0.;
0207 if (f4D) {
0208 for (const auto& tk : *iclus) {
0209 const double time = tk.timeExt();
0210 const double err = tk.dtErrorExt();
0211 const double inverr = err > 0. ? 1.0 / err : 0.;
0212 const double w = inverr * inverr;
0213 sumwt += w * time;
0214 sumwt2 += w * time * time;
0215 sumw += w;
0216 }
0217 meantime = sumwt / sumw;
0218 double sumsq = sumwt2 - sumwt * sumwt / sumw;
0219 double chisq = iclus->size() > 1 ? sumsq / double(iclus->size() - 1) : sumsq / double(iclus->size());
0220 vartime = chisq / sumw;
0221 }
0222
0223 TransientVertex v;
0224 if (algorithm->fitter) {
0225 if (algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) {
0226 v = algorithm->fitter->vertex(*iclus, beamSpot);
0227 } else if (!(algorithm->useBeamConstraint) && (iclus->size() > 1)) {
0228 v = algorithm->fitter->vertex(*iclus);
0229 }
0230 } else if (weightFit) {
0231 std::vector<std::pair<GlobalPoint, GlobalPoint>> points;
0232 if (algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) {
0233 for (const auto& itrack : *iclus) {
0234 GlobalPoint p = itrack.stateAtBeamLine().trackStateAtPCA().position();
0235 GlobalPoint err(itrack.stateAtBeamLine().transverseImpactParameter().error(),
0236 itrack.stateAtBeamLine().transverseImpactParameter().error(),
0237 itrack.track().dzError());
0238 std::pair<GlobalPoint, GlobalPoint> p2(p, err);
0239 points.push_back(p2);
0240 }
0241
0242 v = WeightedMeanFitter::weightedMeanOutlierRejectionBeamSpot(points, *iclus, beamSpot);
0243 if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError))
0244 pvs.push_back(v);
0245 } else if (!(algorithm->useBeamConstraint) && (iclus->size() > 1)) {
0246 for (const auto& itrack : *iclus) {
0247 GlobalPoint p = itrack.impactPointState().globalPosition();
0248 GlobalPoint err(itrack.track().dxyError(), itrack.track().dxyError(), itrack.track().dzError());
0249 std::pair<GlobalPoint, GlobalPoint> p2(p, err);
0250 points.push_back(p2);
0251 }
0252
0253 v = WeightedMeanFitter::weightedMeanOutlierRejection(points, *iclus);
0254 if ((v.positionError().matrix())(2, 2) != (WeightedMeanFitter::startError * WeightedMeanFitter::startError))
0255 pvs.push_back(v);
0256 }
0257 } else
0258 throw VertexException(
0259 "PrimaryVertexProducer: Something went wrong. You are not using the weighted mean fit and no algorithm was "
0260 "selected.");
0261
0262
0263 if (f4D and v.isValid()) {
0264 auto err = v.positionError().matrix4D();
0265 err(3, 3) = vartime;
0266 auto trkWeightMap3d = v.weightMap();
0267 v = TransientVertex(v.position(), meantime, err, v.originalTracks(), v.totalChiSquared(), v.degreesOfFreedom());
0268 v.weightMap(trkWeightMap3d);
0269 }
0270
0271 if (fVerbose) {
0272 if (v.isValid()) {
0273 std::cout << "x,y,z";
0274 if (f4D)
0275 std::cout << ",t";
0276 std::cout << "=" << v.position().x() << " " << v.position().y() << " " << v.position().z();
0277 if (f4D)
0278 std::cout << " " << v.time();
0279 std::cout << " cluster size = " << (*iclus).size() << std::endl;
0280 } else {
0281 std::cout << "Invalid fitted vertex, cluster size=" << (*iclus).size() << std::endl;
0282 }
0283 }
0284
0285
0286 if (v.isValid() && not weightFit && (v.degreesOfFreedom() >= algorithm->minNdof) &&
0287 (!validBS || (*(algorithm->vertexSelector))(v, beamVertexState)))
0288 pvs.push_back(v);
0289 }
0290
0291 if (fVerbose) {
0292 std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << pvs.size() << std::endl;
0293 }
0294
0295 if (clusters.size() > 2 && clusters.size() > 2 * pvs.size())
0296 edm::LogWarning("PrimaryVertexProducer")
0297 << "more than half of candidate vertices lost " << pvs.size() << ' ' << clusters.size();
0298
0299 if (pvs.empty() && seltks.size() > 5)
0300 edm::LogWarning("PrimaryVertexProducer")
0301 << "no vertex found with " << seltks.size() << " tracks and " << clusters.size() << " vertex-candidates";
0302
0303
0304 if (pvs.size() > 1) {
0305 sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
0306 }
0307
0308
0309 for (std::vector<TransientVertex>::const_iterator iv = pvs.begin(); iv != pvs.end(); iv++) {
0310 reco::Vertex v = *iv;
0311 vColl.push_back(v);
0312 }
0313
0314 if (vColl.empty()) {
0315 GlobalError bse(beamSpot.rotatedCovariance3D());
0316 if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
0317 AlgebraicSymMatrix33 we;
0318 we(0, 0) = 10000;
0319 we(1, 1) = 10000;
0320 we(2, 2) = 10000;
0321 vColl.push_back(reco::Vertex(beamSpot.position(), we, 0., 0., 0));
0322 if (fVerbose) {
0323 std::cout << "RecoVertex/PrimaryVertexProducer: "
0324 << "Beamspot with invalid errors " << bse.matrix() << std::endl;
0325 std::cout << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n";
0326 }
0327 } else {
0328 vColl.push_back(reco::Vertex(beamSpot.position(), beamSpot.rotatedCovariance3D(), 0., 0., 0));
0329 if (fVerbose) {
0330 std::cout << "RecoVertex/PrimaryVertexProducer: "
0331 << " will put Vertex derived from BeamSpot into Event.\n";
0332 }
0333 }
0334 }
0335
0336 if (fVerbose) {
0337 int ivtx = 0;
0338 for (reco::VertexCollection::const_iterator v = vColl.begin(); v != vColl.end(); ++v) {
0339 std::cout << "recvtx " << ivtx++ << "#trk " << std::setw(3) << v->tracksSize() << " chi2 " << std::setw(4)
0340 << v->chi2() << " ndof " << std::setw(3) << v->ndof() << " x " << std::setw(6) << v->position().x()
0341 << " dx " << std::setw(6) << v->xError() << " y " << std::setw(6) << v->position().y() << " dy "
0342 << std::setw(6) << v->yError() << " z " << std::setw(6) << v->position().z() << " dz " << std::setw(6)
0343 << v->zError();
0344 if (f4D) {
0345 std::cout << " t " << std::setw(6) << v->t() << " dt " << std::setw(6) << v->tError();
0346 }
0347 std::cout << std::endl;
0348 }
0349 }
0350
0351 iEvent.put(std::move(result), algorithm->label);
0352 }
0353 }
0354
0355 void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0356
0357 edm::ParameterSetDescription desc;
0358 {
0359 edm::ParameterSetDescription vpsd1;
0360 vpsd1.add<double>("maxDistanceToBeam", 1.0);
0361 vpsd1.add<std::string>("algorithm", "AdaptiveVertexFitter");
0362 vpsd1.add<bool>("useBeamConstraint", false);
0363 vpsd1.add<std::string>("label", "");
0364 vpsd1.add<double>("chi2cutoff", 2.5);
0365 vpsd1.add<double>("minNdof", 0.0);
0366 std::vector<edm::ParameterSet> temp1;
0367 temp1.reserve(2);
0368 {
0369 edm::ParameterSet temp2;
0370 temp2.addParameter<double>("maxDistanceToBeam", 1.0);
0371 temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
0372 temp2.addParameter<bool>("useBeamConstraint", false);
0373 temp2.addParameter<std::string>("label", "");
0374 temp2.addParameter<double>("chi2cutoff", 2.5);
0375 temp2.addParameter<double>("minNdof", 0.0);
0376 temp1.push_back(temp2);
0377 }
0378 {
0379 edm::ParameterSet temp2;
0380 temp2.addParameter<double>("maxDistanceToBeam", 1.0);
0381 temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
0382 temp2.addParameter<bool>("useBeamConstraint", true);
0383 temp2.addParameter<std::string>("label", "WithBS");
0384 temp2.addParameter<double>("chi2cutoff", 2.5);
0385 temp2.addParameter<double>("minNdof", 2.0);
0386 temp1.push_back(temp2);
0387 }
0388 desc.addVPSet("vertexCollections", vpsd1, temp1);
0389 }
0390 desc.addUntracked<bool>("verbose", false);
0391 {
0392 edm::ParameterSetDescription psd0;
0393 TrackFilterForPVFinding::fillPSetDescription(psd0);
0394 psd0.add<int>("numTracksThreshold", 0);
0395 psd0.add<int>("maxNumTracksThreshold", 10000000);
0396 psd0.add<double>("minPtTight", 0.0);
0397 desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
0398 }
0399 desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0400 desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
0401 desc.add<edm::InputTag>("TrackTimeResosLabel", edm::InputTag("dummy_default"));
0402 desc.add<edm::InputTag>("TrackTimesLabel", edm::InputTag("dummy_default"));
0403
0404 {
0405 edm::ParameterSetDescription psd0;
0406 {
0407 edm::ParameterSetDescription psd1;
0408 DAClusterizerInZT_vect::fillPSetDescription(psd1);
0409 psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
0410
0411 edm::ParameterSetDescription psd2;
0412 GapClusterizerInZ::fillPSetDescription(psd2);
0413 psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
0414 }
0415 psd0.add<std::string>("algorithm", "DA_vect");
0416 desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
0417 }
0418
0419 desc.add<bool>("isRecoveryIteration", false);
0420 desc.add<edm::InputTag>("recoveryVtxCollection", {""});
0421
0422 descriptions.add("primaryVertexProducer", desc);
0423 }
0424
0425
0426 DEFINE_FWK_MODULE(PrimaryVertexProducer);