PixelVertexFinding/plugins/PixelVertexProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PixelVertexProducer
0004 // Class:      PixelVertexProducer
0005 //
0006 /**\class PixelVertexProducer PixelVertexProducer.h PixelVertexFinding/interface/PixelVertexProducer.h
0007
0008  Description: This produces 1D (z only) primary vertexes using only pixel information.
0009
0010  Implementation:
0011      This producer can use either the Divisive Primary Vertex Finder
0012      or the Histogramming Primary Vertex Finder (currently not
0013      implemented).  It relies on the PixelTripletProducer and
0014      PixelTrackProducer having already been run upstream.   This is
0015      code ported from ORCA originally written by S Cucciarelli, M
0016      Konecki, D Kotlinski.
0017 */
0018 //
0019 // Original Author:  Aaron Dominguez (UNL)
0020 //         Created:  Thu May 25 10:17:32 CDT 2006
0021 //
0022
0023 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0024 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0025 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0027 #include "DataFormats/VertexReco/interface/Vertex.h"
0028 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/EventSetup.h"
0031 #include "FWCore/Framework/interface/Frameworkfwd.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/Framework/interface/stream/EDProducer.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039 #include "RecoVertex/PixelVertexFinding/interface/DivisiveVertexFinder.h"
0040
0041 #include <memory>
0042 #include <string>
0043 #include <cmath>
0044
0045 class PixelVertexProducer : public edm::stream::EDProducer<> {
0046 public:
0047   explicit PixelVertexProducer(const edm::ParameterSet&);
0048   ~PixelVertexProducer() override;
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051   void produce(edm::Event&, const edm::EventSetup&) override;
0052
0053 private:
0054   // ----------member data ---------------------------
0055   // Turn on debug printing if verbose_ > 0
0056   const int verbose_;
0057   // Tracking cuts before sending tracks to vertex algo
0058   const double ptMin_;
0059   const bool method2;
0060   const edm::InputTag trackCollName;
0061   const edm::EDGetTokenT<reco::TrackCollection> token_Tracks;
0062   const edm::EDGetTokenT<reco::BeamSpot> token_BeamSpot;
0063
0064   DivisiveVertexFinder* dvf_;
0065 };
0066
0067 void PixelVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0068   edm::ParameterSetDescription desc;
0069   desc.add<int>("Verbosity", 0);
0070   desc.add<double>("PtMin", 1.0);
0071   desc.add<bool>("Method2", true);
0072   desc.add<edm::InputTag>("TrackCollection", edm::InputTag("pixelTracks"));
0073   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0074   desc.add<std::string>("Finder", "DivisiveVertexFinder");
0075   desc.add<bool>("UseError", true);
0076   desc.add<bool>("WtAverage", true);
0077   desc.add<double>("ZOffset", 5.0);
0078   desc.add<double>("ZSeparation", 0.05);
0079   desc.add<int>("NTrkMin", 2);
0080
0081   {
0082     edm::ParameterSetDescription PVComparerPSet;
0083     PVComparerPSet.add<double>("track_pt_min", 1.0);
0084     PVComparerPSet.add<double>("track_pt_max", 10.0);
0085     PVComparerPSet.add<double>("track_chi2_max", 999999.);
0086     PVComparerPSet.add<double>("track_prob_min", -1.);
0087     desc.addOptional<edm::ParameterSetDescription>("PVcomparer", PVComparerPSet);
0088   }
0089
0090   descriptions.addWithDefaultLabel(desc);
0091 }
0092
0093 PixelVertexProducer::PixelVertexProducer(const edm::ParameterSet& conf)
0094     // 0 silent, 1 chatty, 2 loud
0095     : verbose_(conf.getParameter<int>("Verbosity")),
0096       // 1.0 GeV
0097       ptMin_(conf.getParameter<double>("PtMin")),
0098       method2(conf.getParameter<bool>("Method2")),
0099       trackCollName(conf.getParameter<edm::InputTag>("TrackCollection")),
0100       token_Tracks(consumes<reco::TrackCollection>(trackCollName)),
0101       token_BeamSpot(consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpot"))) {
0102   // Register my product
0103   produces<reco::VertexCollection>();
0104
0105   // Setup shop
0106   std::string finder = conf.getParameter<std::string>("Finder");  // DivisiveVertexFinder
0107   bool useError = conf.getParameter<bool>("UseError");            // true
0108   bool wtAverage = conf.getParameter<bool>("WtAverage");          // true
0109   double zOffset = conf.getParameter<double>("ZOffset");          // 5.0 sigma
0110   double zSeparation = conf.getParameter<double>("ZSeparation");  // 0.05 cm
0111   int ntrkMin = conf.getParameter<int>("NTrkMin");                // 3
0112   // Tracking requirements before sending a track to be considered for vtx
0113
0114   double track_pt_min = ptMin_;
0115   double track_pt_max = 10.;
0116   double track_chi2_max = 9999999.;
0117   double track_prob_min = -1.;
0118
0119   if (conf.exists("PVcomparer")) {
0120     edm::ParameterSet PVcomparerPSet = conf.getParameter<edm::ParameterSet>("PVcomparer");
0121     track_pt_min = PVcomparerPSet.getParameter<double>("track_pt_min");
0122     if (track_pt_min != ptMin_) {
0123       if (track_pt_min < ptMin_)
0124         edm::LogInfo("PixelVertexProducer")
0125             << "minimum track pT setting differs between PixelVertexProducer (" << ptMin_ << ") and PVcomparer ("
0126             << track_pt_min << ") [PVcomparer considers tracks w/ lower threshold than PixelVertexProducer does] !!!";
0127       else
0128         edm::LogInfo("PixelVertexProducer") << "minimum track pT setting differs between PixelVertexProducer ("
0129                                             << ptMin_ << ") and PVcomparer (" << track_pt_min << ") !!!";
0130     }
0131     track_pt_max = PVcomparerPSet.getParameter<double>("track_pt_max");
0132     track_chi2_max = PVcomparerPSet.getParameter<double>("track_chi2_max");
0133     track_prob_min = PVcomparerPSet.getParameter<double>("track_prob_min");
0134   }
0135
0136   if (finder == "DivisiveVertexFinder") {
0137     if (verbose_ > 0)
0138       edm::LogInfo("PixelVertexProducer") << ": Using the DivisiveVertexFinder\n";
0139     dvf_ = new DivisiveVertexFinder(track_pt_min,
0140                                     track_pt_max,
0141                                     track_chi2_max,
0142                                     track_prob_min,
0143                                     zOffset,
0144                                     ntrkMin,
0145                                     useError,
0146                                     zSeparation,
0147                                     wtAverage,
0148                                     verbose_);
0149   } else {  // Finder not supported, or you made a mistake in your request
0150     // throw an exception once I figure out how CMSSW does this
0151   }
0152 }
0153
0154 PixelVertexProducer::~PixelVertexProducer() { delete dvf_; }
0155
0156 void PixelVertexProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0157   // First fish the pixel tracks out of the event
0158   edm::Handle<reco::TrackCollection> trackCollection;
0159   e.getByToken(token_Tracks, trackCollection);
0160   const reco::TrackCollection tracks = *(trackCollection.product());
0161   if (verbose_ > 0)
0162     edm::LogInfo("PixelVertexProducer") << ": Found " << tracks.size() << " tracks in TrackCollection called "
0163                                         << trackCollName << "\n";
0164
0165   // Second, make a collection of pointers to the tracks we want for the vertex finder
0166   reco::TrackRefVector trks;
0167   for (unsigned int i = 0; i < tracks.size(); i++) {
0168     if (tracks[i].pt() > ptMin_)
0169       trks.push_back(reco::TrackRef(trackCollection, i));
0170   }
0171   if (verbose_ > 0)
0172     edm::LogInfo("PixelVertexProducer") << ": Selected " << trks.size() << " of these tracks for vertexing\n";
0173
0174   edm::Handle<reco::BeamSpot> bsHandle;
0175   e.getByToken(token_BeamSpot, bsHandle);
0176   math::XYZPoint myPoint(0., 0., 0.);
0177   if (bsHandle.isValid())
0178     //FIXME: fix last coordinate with vertex.z() at same time
0179     myPoint = math::XYZPoint(bsHandle->x0(), bsHandle->y0(), 0.);
0180
0181   // Third, ship these tracks off to be vertexed
0182   auto vertexes = std::make_unique<reco::VertexCollection>();
0183   bool ok;
0184   if (method2) {
0185     ok = dvf_->findVertexesAlt(trks,  // input
0186                                *vertexes,
0187                                myPoint);  // output
0188     if (verbose_ > 0)
0189       edm::LogInfo("PixelVertexProducer") << "Method2 returned status of " << ok;
0190   } else {
0191     ok = dvf_->findVertexes(trks,        // input
0192                             *vertexes);  // output
0193     if (verbose_ > 0)
0194       edm::LogInfo("PixelVertexProducer") << "Method1 returned status of " << ok;
0195   }
0196
0197   if (verbose_ > 0) {
0198     edm::LogInfo("PixelVertexProducer") << ": Found " << vertexes->size() << " vertexes\n";
0199     for (unsigned int i = 0; i < vertexes->size(); ++i) {
0200       edm::LogInfo("PixelVertexProducer")
0201           << "Vertex number " << i << " has " << (*vertexes)[i].tracksSize() << " tracks with a position of "
0202           << (*vertexes)[i].z() << " +- " << std::sqrt((*vertexes)[i].covariance(2, 2));
0203     }
0204   }
0205
0206   if (bsHandle.isValid()) {
0207     const reco::BeamSpot& bs = *bsHandle;
0208
0209     for (unsigned int i = 0; i < vertexes->size(); ++i) {
0210       double z = (*vertexes)[i].z();
0211       double x = bs.x0() + bs.dxdz() * (z - bs.z0());
0212       double y = bs.y0() + bs.dydz() * (z - bs.z0());
0213       reco::Vertex v(reco::Vertex::Point(x, y, z),
0214                      (*vertexes)[i].error(),
0215                      (*vertexes)[i].chi2(),
0216                      (*vertexes)[i].ndof(),
0217                      (*vertexes)[i].tracksSize());
0218       //Copy also the tracks
0219       for (std::vector<reco::TrackBaseRef>::const_iterator it = (*vertexes)[i].tracks_begin();
0220            it != (*vertexes)[i].tracks_end();
0221            it++) {
0222         v.add(*it);
0223       }
0224       (*vertexes)[i] = v;
0225     }
0226   } else {
0227     edm::LogWarning("PixelVertexProducer") << "No beamspot found. Using returning vertexes with (0,0,Z) ";
0228   }
0229
0230   if (vertexes->empty() && bsHandle.isValid()) {
0231     const reco::BeamSpot& bs = *bsHandle;
0232
0233     GlobalError bse(bs.rotatedCovariance3D());
0234     if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
0235       AlgebraicSymMatrix33 we;
0236       we(0, 0) = 10000;
0237       we(1, 1) = 10000;
0238       we(2, 2) = 10000;
0239       vertexes->push_back(reco::Vertex(bs.position(), we, 0., 0., 0));
0240
0241       edm::LogInfo("PixelVertexProducer")
0242           << "No vertices found. Beamspot with invalid errors " << bse.matrix() << std::endl
0243           << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n"
0244           << (*vertexes)[0].x() << "\n"
0245           << (*vertexes)[0].y() << "\n"
0246           << (*vertexes)[0].z() << "\n";
0247     } else {
0248       vertexes->push_back(reco::Vertex(bs.position(), bs.rotatedCovariance3D(), 0., 0., 0));
0249
0250       edm::LogInfo("PixelVertexProducer") << "No vertices found. Will put Vertex derived from BeamSpot into Event:\n"
0251                                           << (*vertexes)[0].x() << "\n"
0252                                           << (*vertexes)[0].y() << "\n"
0253                                           << (*vertexes)[0].z() << "\n";
0254     }
0255   }
0256
0257   else if (vertexes->empty() && !bsHandle.isValid()) {
0258     edm::LogWarning("PixelVertexProducer") << "No beamspot and no vertex found. No vertex returned.";
0259   }
0260
0261   e.put(std::move(vertexes));
0262 }
0263
0264 DEFINE_FWK_MODULE(PixelVertexProducer);