Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // select and configure the track selection
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   // select and configure the track clusterizer
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   // provide the vectorized version of the clusterizer, if supported by the build
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   // select and configure the vertex fitters
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   //check if this is a recovery iteration
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   // get the BeamSpot, it will always be needed, even when not used as a constraint
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   //if this is a recovery iteration, check if we already have a valid PV
0143   if (fRecoveryIteration) {
0144     auto const& oldVertices = iEvent.get(recoveryVtxToken);
0145     //look for the first valid (not-BeamSpot) vertex
0146     for (auto const& old : oldVertices) {
0147       if (!(old.isFake())) {
0148         //found a valid vertex, write the first one to the collection and return
0149         //otherwise continue with regular vertexing procedure
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   // get RECO tracks from the event
0159   // `tks` can be used as a ptr to a reco::TrackCollection
0160   edm::Handle<reco::TrackCollection> tks;
0161   iEvent.getByToken(trkToken, tks);
0162 
0163   // interface RECO tracks to vertex reconstruction
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   // select tracks
0183   std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
0184 
0185   // clusterize tracks in Z
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   // vertex fits
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         }  // else: no fit ==> v.isValid()=False
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);  //FIX with constants
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       // 4D vertices: add timing information
0263       if (f4D and v.isValid()) {
0264         auto err = v.positionError().matrix4D();
0265         err(3, 3) = vartime;
0266         auto trkWeightMap3d = v.weightMap();  // copy the 3d-fit weights
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       //for weightFit we have already pushed it above (no timing infomration anyway)
0286       if (v.isValid() && not weightFit && (v.degreesOfFreedom() >= algorithm->minNdof) &&
0287           (!validBS || (*(algorithm->vertexSelector))(v, beamVertexState)))
0288         pvs.push_back(v);
0289     }  // end of cluster loop
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     // sort vertices by pt**2  vertex (aka signal vertex tagging)
0304     if (pvs.size() > 1) {
0305       sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
0306     }
0307 
0308     // convert transient vertices returned by the theAlgo to (reco) vertices
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   // offlinePrimaryVertices
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);            // HI only
0395     psd0.add<int>("maxNumTracksThreshold", 10000000);  // HI only
0396     psd0.add<double>("minPtTight", 0.0);               // HI only
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"));  // 4D only
0402   desc.add<edm::InputTag>("TrackTimesLabel", edm::InputTag("dummy_default"));      // 4D only
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 //define this as a plug-in
0426 DEFINE_FWK_MODULE(PrimaryVertexProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!