File indexing completed on 2023-10-25 10:03:27
0001
0002 #include "RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducerAlgorithm.h"
0003
0004 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0005 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0011 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0012 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0013
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0016 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0017 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0018
0019 PrimaryVertexProducerAlgorithm::PrimaryVertexProducerAlgorithm(const edm::ParameterSet& conf) : theConfig(conf) {
0020 fVerbose = conf.getUntrackedParameter<bool>("verbose", false);
0021 trackLabel = conf.getParameter<edm::InputTag>("TrackLabel");
0022 beamSpotLabel = conf.getParameter<edm::InputTag>("beamSpotLabel");
0023
0024
0025 std::string trackSelectionAlgorithm =
0026 conf.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
0027 if (trackSelectionAlgorithm == "filter") {
0028 theTrackFilter = new TrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0029 } else if (trackSelectionAlgorithm == "filterWithThreshold") {
0030 theTrackFilter = new HITrackFilterForPVFinding(conf.getParameter<edm::ParameterSet>("TkFilterParameters"));
0031 } else {
0032 throw VertexException("PrimaryVertexProducerAlgorithm: unknown track selection algorithm: " +
0033 trackSelectionAlgorithm);
0034 }
0035
0036
0037 std::string clusteringAlgorithm =
0038 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
0039 if (clusteringAlgorithm == "gap") {
0040 theTrackClusterizer = new GapClusterizerInZ(
0041 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
0042 } else if (clusteringAlgorithm == "DA") {
0043 theTrackClusterizer = new DAClusterizerInZ(
0044 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0045 }
0046
0047 else if (clusteringAlgorithm == "DA_vect") {
0048 theTrackClusterizer = new DAClusterizerInZ_vect(
0049 conf.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
0050 }
0051
0052 else {
0053 throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
0054 }
0055
0056
0057 std::vector<edm::ParameterSet> vertexCollections =
0058 conf.getParameter<std::vector<edm::ParameterSet> >("vertexCollections");
0059
0060 for (std::vector<edm::ParameterSet>::const_iterator algoconf = vertexCollections.begin();
0061 algoconf != vertexCollections.end();
0062 algoconf++) {
0063 algo algorithm;
0064 std::string fitterAlgorithm = algoconf->getParameter<std::string>("algorithm");
0065 if (fitterAlgorithm == "KalmanVertexFitter") {
0066 algorithm.fitter = new KalmanVertexFitter();
0067 } else if (fitterAlgorithm == "AdaptiveVertexFitter") {
0068 algorithm.fitter = new AdaptiveVertexFitter();
0069 } else {
0070 throw VertexException("PrimaryVertexProducerAlgorithm: unknown algorithm: " + fitterAlgorithm);
0071 }
0072 algorithm.label = algoconf->getParameter<std::string>("label");
0073 algorithm.minNdof = algoconf->getParameter<double>("minNdof");
0074 algorithm.useBeamConstraint = algoconf->getParameter<bool>("useBeamConstraint");
0075 algorithm.vertexSelector =
0076 new VertexCompatibleWithBeam(VertexDistanceXY(), algoconf->getParameter<double>("maxDistanceToBeam"));
0077 algorithms.push_back(algorithm);
0078 }
0079 }
0080
0081 PrimaryVertexProducerAlgorithm::~PrimaryVertexProducerAlgorithm() {
0082 if (theTrackFilter)
0083 delete theTrackFilter;
0084 if (theTrackClusterizer)
0085 delete theTrackClusterizer;
0086 for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0087 if (algorithm->fitter)
0088 delete algorithm->fitter;
0089 if (algorithm->vertexSelector)
0090 delete algorithm->vertexSelector;
0091 }
0092 }
0093
0094
0095
0096
0097
0098
0099 std::vector<TransientVertex> PrimaryVertexProducerAlgorithm::vertices(
0100 const std::vector<reco::TransientTrack>& tracks) const {
0101 throw VertexException("PrimaryVertexProducerAlgorithm: cannot make a Primary Vertex without a beam spot");
0102
0103 return std::vector<TransientVertex>();
0104 }
0105
0106 std::vector<TransientVertex> PrimaryVertexProducerAlgorithm::vertices(const std::vector<reco::TransientTrack>& t_tks,
0107 const reco::BeamSpot& beamSpot,
0108 const std::string& label) const {
0109 bool validBS = true;
0110 VertexState beamVertexState(beamSpot);
0111 if ((beamVertexState.error().cxx() <= 0.) || (beamVertexState.error().cyy() <= 0.) ||
0112 (beamVertexState.error().czz() <= 0.)) {
0113 validBS = false;
0114 edm::LogError("UnusableBeamSpot") << "Beamspot with invalid errors " << beamVertexState.error().matrix();
0115 }
0116
0117
0118
0119
0120
0121
0122
0123 std::vector<reco::TransientTrack> seltks = theTrackFilter->select(t_tks);
0124
0125
0126 std::vector<std::vector<reco::TransientTrack> > clusters = theTrackClusterizer->clusterize(seltks);
0127 if (fVerbose) {
0128 std::cout << " clustering returned " << clusters.size() << " clusters from " << seltks.size()
0129 << " selected tracks" << std::endl;
0130 }
0131
0132
0133 for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0134 if (!(algorithm->label == label))
0135 continue;
0136
0137
0138
0139
0140 std::vector<TransientVertex> pvs;
0141 for (std::vector<std::vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin();
0142 iclus != clusters.end();
0143 iclus++) {
0144 TransientVertex v;
0145 if (algorithm->useBeamConstraint && validBS && ((*iclus).size() > 1)) {
0146 v = algorithm->fitter->vertex(*iclus, beamSpot);
0147
0148 } else if (!(algorithm->useBeamConstraint) && ((*iclus).size() > 1)) {
0149 v = algorithm->fitter->vertex(*iclus);
0150
0151 }
0152
0153 if (fVerbose) {
0154 if (v.isValid())
0155 std::cout << "x,y,z=" << v.position().x() << " " << v.position().y() << " " << v.position().z() << std::endl;
0156 else
0157 std::cout << "Invalid fitted vertex\n";
0158 }
0159
0160 if (v.isValid() && (v.degreesOfFreedom() >= algorithm->minNdof) &&
0161 (!validBS || (*(algorithm->vertexSelector))(v, beamVertexState)))
0162 pvs.push_back(v);
0163 }
0164
0165 if (fVerbose) {
0166 std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << pvs.size() << std::endl;
0167 }
0168
0169
0170 if (pvs.size() > 1) {
0171 sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
0172 }
0173
0174 return pvs;
0175 }
0176
0177 std::vector<TransientVertex> dummy;
0178 return dummy;
0179 }