Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-23 22:40:21

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 
0026   // select and configure the track selection
0027   std::string trackSelectionAlgorithm =
0028       conf.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
0029   if (trackSelectionAlgorithm == "filter") {
0030     theTrackFilter = new TrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0031   } else if (trackSelectionAlgorithm == "filterWithThreshold") {
0032     theTrackFilter = new HITrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0033   } else {
0034     throw VertexException("PrimaryVertexProducer: unknown track selection algorithm: " + trackSelectionAlgorithm);
0035   }
0036 
0037   // select and configure the track clusterizer
0038   std::string clusteringAlgorithm =
0039       conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
0040   if (clusteringAlgorithm == "gap") {
0041     theTrackClusterizer = new GapClusterizerInZ(
0042         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
0043   } else if (clusteringAlgorithm == "DA") {
0044     theTrackClusterizer = new DAClusterizerInZ(
0045         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0046   }
0047   // provide the vectorized version of the clusterizer, if supported by the build
0048   else if (clusteringAlgorithm == "DA_vect") {
0049     theTrackClusterizer = new DAClusterizerInZ_vect(
0050         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0051   } else if (clusteringAlgorithm == "DA2D_vect") {
0052     theTrackClusterizer = new DAClusterizerInZT_vect(
0053         conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0054     f4D = true;
0055   }
0056 
0057   else {
0058     throw VertexException("PrimaryVertexProducer: unknown clustering algorithm: " + clusteringAlgorithm);
0059   }
0060 
0061   if (f4D) {
0062     trkTimesToken = consumes<edm::ValueMap<float> >(conf.getParameter<edm::InputTag>("TrackTimesLabel"));
0063     trkTimeResosToken = consumes<edm::ValueMap<float> >(conf.getParameter<edm::InputTag>("TrackTimeResosLabel"));
0064   }
0065 
0066   // select and configure the vertex fitters
0067   std::vector<edm::ParameterSet> vertexCollections =
0068       conf.getParameter<std::vector<edm::ParameterSet> >("vertexCollections");
0069 
0070   for (std::vector<edm::ParameterSet>::const_iterator algoconf = vertexCollections.begin();
0071        algoconf != vertexCollections.end();
0072        algoconf++) {
0073     algo algorithm;
0074     std::string fitterAlgorithm = algoconf->getParameter<std::string>("algorithm");
0075     if (fitterAlgorithm == "KalmanVertexFitter") {
0076       algorithm.fitter = new KalmanVertexFitter();
0077     } else if (fitterAlgorithm == "AdaptiveVertexFitter") {
0078       algorithm.fitter = new AdaptiveVertexFitter(GeometricAnnealing(algoconf->getParameter<double>("chi2cutoff")));
0079     } else {
0080       throw VertexException("PrimaryVertexProducer: unknown algorithm: " + fitterAlgorithm);
0081     }
0082     algorithm.label = algoconf->getParameter<std::string>("label");
0083     algorithm.minNdof = algoconf->getParameter<double>("minNdof");
0084     algorithm.useBeamConstraint = algoconf->getParameter<bool>("useBeamConstraint");
0085     algorithm.vertexSelector =
0086         new VertexCompatibleWithBeam(VertexDistanceXY(), algoconf->getParameter<double>("maxDistanceToBeam"));
0087     algorithms.push_back(algorithm);
0088 
0089     produces<reco::VertexCollection>(algorithm.label);
0090   }
0091 
0092   //check if this is a recovery iteration
0093   fRecoveryIteration = conf.getParameter<bool>("isRecoveryIteration");
0094   if (fRecoveryIteration) {
0095     if (algorithms.empty()) {
0096       throw VertexException("PrimaryVertexProducer: No algorithm specified. ");
0097     } else if (algorithms.size() > 1) {
0098       throw VertexException(
0099           "PrimaryVertexProducer: Running in Recovery mode and more than one algorithm specified.  Please "
0100           "only one algorithm.");
0101     }
0102     recoveryVtxToken = consumes<reco::VertexCollection>(conf.getParameter<edm::InputTag>("recoveryVtxCollection"));
0103   }
0104 }
0105 
0106 PrimaryVertexProducer::~PrimaryVertexProducer() {
0107   if (theTrackFilter)
0108     delete theTrackFilter;
0109   if (theTrackClusterizer)
0110     delete theTrackClusterizer;
0111   for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0112     if (algorithm->fitter)
0113       delete algorithm->fitter;
0114     if (algorithm->vertexSelector)
0115       delete algorithm->vertexSelector;
0116   }
0117 }
0118 
0119 void PrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0120   // get the BeamSpot, it will always be needed, even when not used as a constraint
0121   reco::BeamSpot beamSpot;
0122   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0123   iEvent.getByToken(bsToken, recoBeamSpotHandle);
0124   if (recoBeamSpotHandle.isValid()) {
0125     beamSpot = *recoBeamSpotHandle;
0126   } else {
0127     edm::LogError("UnusableBeamSpot") << "No beam spot available from EventSetup";
0128   }
0129 
0130   bool validBS = true;
0131   VertexState beamVertexState(beamSpot);
0132   if ((beamVertexState.error().cxx() <= 0.) || (beamVertexState.error().cyy() <= 0.) ||
0133       (beamVertexState.error().czz() <= 0.)) {
0134     validBS = false;
0135     edm::LogError("UnusableBeamSpot") << "Beamspot with invalid errors " << beamVertexState.error().matrix();
0136   }
0137 
0138   //if this is a recovery iteration, check if we already have a valid PV
0139   if (fRecoveryIteration) {
0140     auto const& oldVertices = iEvent.get(recoveryVtxToken);
0141     //look for the first valid (not-BeamSpot) vertex
0142     for (auto const& old : oldVertices) {
0143       if (!(old.isFake())) {
0144         //found a valid vertex, write the first one to the collection and return
0145         //otherwise continue with regular vertexing procedure
0146         auto result = std::make_unique<reco::VertexCollection>();
0147         result->push_back(old);
0148         iEvent.put(std::move(result), algorithms.begin()->label);
0149         return;
0150       }
0151     }
0152   }
0153 
0154   // get RECO tracks from the event
0155   // `tks` can be used as a ptr to a reco::TrackCollection
0156   edm::Handle<reco::TrackCollection> tks;
0157   iEvent.getByToken(trkToken, tks);
0158 
0159   // interface RECO tracks to vertex reconstruction
0160   const auto& theB = &iSetup.getData(theTTBToken);
0161   std::vector<reco::TransientTrack> t_tks;
0162 
0163   if (f4D) {
0164     edm::Handle<edm::ValueMap<float> > trackTimesH;
0165     edm::Handle<edm::ValueMap<float> > trackTimeResosH;
0166     iEvent.getByToken(trkTimesToken, trackTimesH);
0167     iEvent.getByToken(trkTimeResosToken, trackTimeResosH);
0168     t_tks = (*theB).build(tks, beamSpot, *(trackTimesH.product()), *(trackTimeResosH.product()));
0169   } else {
0170     t_tks = (*theB).build(tks, beamSpot);
0171   }
0172   if (fVerbose) {
0173     std::cout << "RecoVertex/PrimaryVertexProducer"
0174               << "Found: " << t_tks.size() << " reconstructed tracks"
0175               << "\n";
0176   }
0177 
0178   // select tracks
0179   std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
0180 
0181   // clusterize tracks in Z
0182   std::vector<std::vector<reco::TransientTrack> >&& clusters = theTrackClusterizer->clusterize(seltks);
0183 
0184   if (fVerbose) {
0185     std::cout << " clustering returned  " << clusters.size() << " clusters  from " << seltks.size()
0186               << " selected tracks" << std::endl;
0187   }
0188 
0189   // vertex fits
0190   for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0191     auto result = std::make_unique<reco::VertexCollection>();
0192     reco::VertexCollection& vColl = (*result);
0193 
0194     std::vector<TransientVertex> pvs;
0195     for (std::vector<std::vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin();
0196          iclus != clusters.end();
0197          iclus++) {
0198       double sumwt = 0.;
0199       double sumwt2 = 0.;
0200       double sumw = 0.;
0201       double meantime = 0.;
0202       double vartime = 0.;
0203       if (f4D) {
0204         for (const auto& tk : *iclus) {
0205           const double time = tk.timeExt();
0206           const double err = tk.dtErrorExt();
0207           const double inverr = err > 0. ? 1.0 / err : 0.;
0208           const double w = inverr * inverr;
0209           sumwt += w * time;
0210           sumwt2 += w * time * time;
0211           sumw += w;
0212         }
0213         meantime = sumwt / sumw;
0214         double sumsq = sumwt2 - sumwt * sumwt / sumw;
0215         double chisq = iclus->size() > 1 ? sumsq / double(iclus->size() - 1) : sumsq / double(iclus->size());
0216         vartime = chisq / sumw;
0217       }
0218 
0219       TransientVertex v;
0220       if (algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) {
0221         v = algorithm->fitter->vertex(*iclus, beamSpot);
0222       } else if (!(algorithm->useBeamConstraint) && (iclus->size() > 1)) {
0223         v = algorithm->fitter->vertex(*iclus);
0224       }  // else: no fit ==> v.isValid()=False
0225 
0226       // 4D vertices: add timing information
0227       if (f4D and v.isValid()) {
0228         auto err = v.positionError().matrix4D();
0229         err(3, 3) = vartime;
0230         auto trkWeightMap3d = v.weightMap();  // copy the 3d-fit weights
0231         v = TransientVertex(v.position(), meantime, err, v.originalTracks(), v.totalChiSquared(), v.degreesOfFreedom());
0232         v.weightMap(trkWeightMap3d);
0233       }
0234 
0235       if (fVerbose) {
0236         if (v.isValid()) {
0237           std::cout << "x,y,z";
0238           if (f4D)
0239             std::cout << ",t";
0240           std::cout << "=" << v.position().x() << " " << v.position().y() << " " << v.position().z();
0241           if (f4D)
0242             std::cout << " " << v.time();
0243           std::cout << " cluster size = " << (*iclus).size() << std::endl;
0244         } else {
0245           std::cout << "Invalid fitted vertex,  cluster size=" << (*iclus).size() << std::endl;
0246         }
0247       }
0248 
0249       if (v.isValid() && (v.degreesOfFreedom() >= algorithm->minNdof) &&
0250           (!validBS || (*(algorithm->vertexSelector))(v, beamVertexState)))
0251         pvs.push_back(v);
0252     }  // end of cluster loop
0253 
0254     if (fVerbose) {
0255       std::cout << "PrimaryVertexProducerAlgorithm::vertices  candidates =" << pvs.size() << std::endl;
0256     }
0257 
0258     if (clusters.size() > 2 && clusters.size() > 2 * pvs.size())
0259       edm::LogWarning("PrimaryVertexProducer")
0260           << "more than half of candidate vertices lost " << pvs.size() << ' ' << clusters.size();
0261 
0262     if (pvs.empty() && seltks.size() > 5)
0263       edm::LogWarning("PrimaryVertexProducer")
0264           << "no vertex found with " << seltks.size() << " tracks and " << clusters.size() << " vertex-candidates";
0265 
0266     // sort vertices by pt**2  vertex (aka signal vertex tagging)
0267     if (pvs.size() > 1) {
0268       sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
0269     }
0270 
0271     // convert transient vertices returned by the theAlgo to (reco) vertices
0272     for (std::vector<TransientVertex>::const_iterator iv = pvs.begin(); iv != pvs.end(); iv++) {
0273       reco::Vertex v = *iv;
0274       vColl.push_back(v);
0275     }
0276 
0277     if (vColl.empty()) {
0278       GlobalError bse(beamSpot.rotatedCovariance3D());
0279       if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
0280         AlgebraicSymMatrix33 we;
0281         we(0, 0) = 10000;
0282         we(1, 1) = 10000;
0283         we(2, 2) = 10000;
0284         vColl.push_back(reco::Vertex(beamSpot.position(), we, 0., 0., 0));
0285         if (fVerbose) {
0286           std::cout << "RecoVertex/PrimaryVertexProducer: "
0287                     << "Beamspot with invalid errors " << bse.matrix() << std::endl;
0288           std::cout << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n";
0289         }
0290       } else {
0291         vColl.push_back(reco::Vertex(beamSpot.position(), beamSpot.rotatedCovariance3D(), 0., 0., 0));
0292         if (fVerbose) {
0293           std::cout << "RecoVertex/PrimaryVertexProducer: "
0294                     << " will put Vertex derived from BeamSpot into Event.\n";
0295         }
0296       }
0297     }
0298 
0299     if (fVerbose) {
0300       int ivtx = 0;
0301       for (reco::VertexCollection::const_iterator v = vColl.begin(); v != vColl.end(); ++v) {
0302         std::cout << "recvtx " << ivtx++ << "#trk " << std::setw(3) << v->tracksSize() << " chi2 " << std::setw(4)
0303                   << v->chi2() << " ndof " << std::setw(3) << v->ndof() << " x " << std::setw(6) << v->position().x()
0304                   << " dx " << std::setw(6) << v->xError() << " y " << std::setw(6) << v->position().y() << " dy "
0305                   << std::setw(6) << v->yError() << " z " << std::setw(6) << v->position().z() << " dz " << std::setw(6)
0306                   << v->zError();
0307         if (f4D) {
0308           std::cout << " t " << std::setw(6) << v->t() << " dt " << std::setw(6) << v->tError();
0309         }
0310         std::cout << std::endl;
0311       }
0312     }
0313 
0314     iEvent.put(std::move(result), algorithm->label);
0315   }
0316 }
0317 
0318 void PrimaryVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0319   // offlinePrimaryVertices
0320   edm::ParameterSetDescription desc;
0321   {
0322     edm::ParameterSetDescription vpsd1;
0323     vpsd1.add<double>("maxDistanceToBeam", 1.0);
0324     vpsd1.add<std::string>("algorithm", "AdaptiveVertexFitter");
0325     vpsd1.add<bool>("useBeamConstraint", false);
0326     vpsd1.add<std::string>("label", "");
0327     vpsd1.add<double>("chi2cutoff", 2.5);
0328     vpsd1.add<double>("minNdof", 0.0);
0329     std::vector<edm::ParameterSet> temp1;
0330     temp1.reserve(2);
0331     {
0332       edm::ParameterSet temp2;
0333       temp2.addParameter<double>("maxDistanceToBeam", 1.0);
0334       temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
0335       temp2.addParameter<bool>("useBeamConstraint", false);
0336       temp2.addParameter<std::string>("label", "");
0337       temp2.addParameter<double>("chi2cutoff", 2.5);
0338       temp2.addParameter<double>("minNdof", 0.0);
0339       temp1.push_back(temp2);
0340     }
0341     {
0342       edm::ParameterSet temp2;
0343       temp2.addParameter<double>("maxDistanceToBeam", 1.0);
0344       temp2.addParameter<std::string>("algorithm", "AdaptiveVertexFitter");
0345       temp2.addParameter<bool>("useBeamConstraint", true);
0346       temp2.addParameter<std::string>("label", "WithBS");
0347       temp2.addParameter<double>("chi2cutoff", 2.5);
0348       temp2.addParameter<double>("minNdof", 2.0);
0349       temp1.push_back(temp2);
0350     }
0351     desc.addVPSet("vertexCollections", vpsd1, temp1);
0352   }
0353   desc.addUntracked<bool>("verbose", false);
0354   {
0355     edm::ParameterSetDescription psd0;
0356     TrackFilterForPVFinding::fillPSetDescription(psd0);
0357     psd0.add<int>("numTracksThreshold", 0);            // HI only
0358     psd0.add<int>("maxNumTracksThreshold", 10000000);  // HI only
0359     psd0.add<double>("minPtTight", 0.0);               // HI only
0360     desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
0361   }
0362   desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0363   desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
0364   desc.add<edm::InputTag>("TrackTimeResosLabel", edm::InputTag("dummy_default"));  // 4D only
0365   desc.add<edm::InputTag>("TrackTimesLabel", edm::InputTag("dummy_default"));      // 4D only
0366 
0367   {
0368     edm::ParameterSetDescription psd0;
0369     {
0370       edm::ParameterSetDescription psd1;
0371       DAClusterizerInZT_vect::fillPSetDescription(psd1);
0372       psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
0373 
0374       edm::ParameterSetDescription psd2;
0375       GapClusterizerInZ::fillPSetDescription(psd2);
0376       psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
0377     }
0378     psd0.add<std::string>("algorithm", "DA_vect");
0379     desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
0380   }
0381 
0382   desc.add<bool>("isRecoveryIteration", false);
0383   desc.add<edm::InputTag>("recoveryVtxCollection", {""});
0384 
0385   descriptions.add("primaryVertexProducer", desc);
0386 }
0387 
0388 //define this as a plug-in
0389 DEFINE_FWK_MODULE(PrimaryVertexProducer);