PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc

0001 ///////////////   OBSOLETE ////////////////////
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   // select and configure the track selection
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   // select and configure the track clusterizer
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   // provide the vectorized version of the clusterizer, if supported by the build
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   // select and configure the vertex fitters
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 // member functions
0096 //
0097
0098 // obsolete method, unfortunately required through inheritance from  VertexReconstructor
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   //   // get RECO tracks from the event
0118   //   // `tks` can be used as a ptr to a reco::TrackCollection
0119   //   edm::Handle<reco::TrackCollection> tks;
0120   //   iEvent.getByLabel(trackLabel, tks);
0121
0122   // select tracks
0123   std::vector<reco::TransientTrack> seltks = theTrackFilter->select(t_tks);
0124
0125   // clusterize tracks in Z
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   // vertex fits
0133   for (std::vector<algo>::const_iterator algorithm = algorithms.begin(); algorithm != algorithms.end(); algorithm++) {
0134     if (!(algorithm->label == label))
0135       continue;
0136
0137     //std::auto_ptr<reco::VertexCollection> result(new reco::VertexCollection);
0138     // reco::VertexCollection vColl;
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       }  // else: no fit ==> v.isValid()=False
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     }  // end of cluster loop
0164
0165     if (fVerbose) {
0166       std::cout << "PrimaryVertexProducerAlgorithm::vertices  candidates =" << pvs.size() << std::endl;
0167     }
0168
0169     // sort vertices by pt**2  vertex (aka signal vertex tagging)
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;  //avoid compiler warning, should never be here
0179 }