PixelVertexFinding/plugins/FastPrimaryVertexProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FastPrimaryVertexProducer
0004 // Class:      FastPrimaryVertexProducer
0005 //
0006 /**\class FastPrimaryVertexProducer FastPrimaryVertexProducer.cc RecoBTag/FastPrimaryVertexProducer/src/FastPrimaryVertexProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Andrea RIZZI
0015 //         Created:  Thu Dec 22 14:51:44 CET 2011
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <vector>
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDProducer.h"
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0031 #include "DataFormats/JetReco/interface/CaloJet.h"
0032 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0033
0034 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0036
0037 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0038 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0039 #include "DataFormats/Math/interface/deltaPhi.h"
0040 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0041 #include "DataFormats/VertexReco/interface/Vertex.h"
0042 #include "DataFormats/TrackReco/interface/Track.h"
0043 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0044 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0045 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0046
0047 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0048 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0049 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0050
0051 #include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
0052 #include "CommonTools/Clustering1D/interface/Cluster1DMerger.h"
0053 #include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
0054 #define HaveMtv
0055 #define HaveFsmw
0056 #define HaveDivisive
0057 #ifdef HaveMtv
0058 #include "CommonTools/Clustering1D/interface/MtvClusterizer1D.h"
0059 #endif
0060 #ifdef HaveFsmw
0061 #include "CommonTools/Clustering1D/interface/FsmwClusterizer1D.h"
0062 #endif
0063 #ifdef HaveDivisive
0064 #include "CommonTools/Clustering1D/interface/DivisiveClusterizer1D.h"
0065 #endif
0066
0067 using namespace std;
0068
0069 //
0070 // class declaration
0071 //
0072
0073 class FastPrimaryVertexProducer : public edm::global::EDProducer<> {
0074 public:
0075   explicit FastPrimaryVertexProducer(const edm::ParameterSet&);
0076
0077 private:
0078   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0079
0080   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> const m_geomToken;
0081   edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> const m_pixelCPEToken;
0082   edm::EDGetTokenT<SiPixelClusterCollectionNew> m_clusters;
0083   edm::EDGetTokenT<edm::View<reco::Jet> > m_jets;
0084   edm::EDGetTokenT<reco::BeamSpot> m_beamSpot;
0085   double m_maxZ;
0086   double m_maxSizeX;
0087   double m_maxDeltaPhi;
0088   double m_clusterLength;
0089 };
0090
0091 FastPrimaryVertexProducer::FastPrimaryVertexProducer(const edm::ParameterSet& iConfig)
0092     : m_geomToken(esConsumes()),
0093       m_pixelCPEToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCPE")))) {
0094   m_clusters = consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusters"));
0095   m_jets = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
0096   m_beamSpot = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"));
0097   m_maxZ = iConfig.getParameter<double>("maxZ");
0098   m_maxSizeX = iConfig.getParameter<double>("maxSizeX");
0099   m_maxDeltaPhi = iConfig.getParameter<double>("maxDeltaPhi");
0100   m_clusterLength = iConfig.getParameter<double>("clusterLength");
0101   produces<reco::VertexCollection>();
0102 }
0103
0104 void FastPrimaryVertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0105   using namespace edm;
0106   using namespace reco;
0107   using namespace std;
0108
0109   Handle<SiPixelClusterCollectionNew> cH;
0110   iEvent.getByToken(m_clusters, cH);
0111   const SiPixelClusterCollectionNew& pixelClusters = *cH.product();
0112
0113   Handle<edm::View<reco::Jet> > jH;
0114   iEvent.getByToken(m_jets, jH);
0115   const edm::View<reco::Jet>& jets = *jH.product();
0116
0117   CaloJetCollection selectedJets;
0118   for (edm::View<reco::Jet>::const_iterator it = jets.begin(); it != jets.end(); it++) {
0119     if (it->pt() > 40 && fabs(it->eta()) < 1.6) {
0120       const CaloJet* ca = dynamic_cast<const CaloJet*>(&(*it));
0121       if (ca == nullptr)
0122         abort();
0123       selectedJets.push_back(*ca);
0124       //    std::cout << "Jet eta,phi,pt: "<< it->eta() << "," << it->phi() << "," << it->pt()   << std::endl;
0125     }
0126   }
0127
0128   const PixelClusterParameterEstimator* pp = &iSetup.getData(m_pixelCPEToken);
0129
0130   edm::Handle<BeamSpot> beamSpot;
0131   iEvent.getByToken(m_beamSpot, beamSpot);
0132
0133   const TrackerGeometry* trackerGeometry = &iSetup.getData(m_geomToken);
0134
0135   float lengthBmodule = 6.66;  //cm
0136   std::vector<float> zProjections;
0137   for (CaloJetCollection::const_iterator jit = selectedJets.begin(); jit != selectedJets.end(); jit++) {
0138     float px = jit->px();
0139     float py = jit->py();
0140     float pz = jit->pz();
0141     float pt = jit->pt();
0142
0143     float jetZOverRho = jit->momentum().Z() / jit->momentum().Rho();
0144     int minSizeY = fabs(2. * jetZOverRho) - 1;
0145     int maxSizeY = fabs(2. * jetZOverRho) + 2;
0146     if (fabs(jit->eta()) > 1.6) {
0147       minSizeY = 1;
0148     }
0149
0150     for (SiPixelClusterCollectionNew::const_iterator it = pixelClusters.begin(); it != pixelClusters.end();
0151          it++)  //Loop on pixel modules with clusters
0152     {
0153       DetId id = it->detId();
0154       const edmNew::DetSet<SiPixelCluster>& detset = (*it);
0155       Point3DBase<float, GlobalTag> modulepos = trackerGeometry->idToDet(id)->position();
0156       float zmodule = modulepos.z() -
0157                       ((modulepos.x() - beamSpot->x0()) * px + (modulepos.y() - beamSpot->y0()) * py) / pt * pz / pt;
0158       if ((fabs(deltaPhi(jit->momentum().Phi(), modulepos.phi())) < m_maxDeltaPhi * 2) &&
0159           (fabs(zmodule) < (m_maxZ + lengthBmodule / 2))) {
0160         for (size_t j = 0; j < detset.size(); j++)  // Loop on pixel clusters on this module
0161         {
0162           const SiPixelCluster& aCluster = detset[j];
0163           if (aCluster.sizeX() < m_maxSizeX && aCluster.sizeY() >= minSizeY && aCluster.sizeY() <= maxSizeY) {
0164             Point3DBase<float, GlobalTag> v = trackerGeometry->idToDet(id)->surface().toGlobal(
0165                 pp->localParametersV(aCluster, (*trackerGeometry->idToDetUnit(id)))[0].first);
0166             GlobalPoint v_bs(v.x() - beamSpot->x0(), v.y() - beamSpot->y0(), v.z());
0167             if (fabs(deltaPhi(jit->momentum().Phi(), v_bs.phi())) < m_maxDeltaPhi) {
0168               float z = v.z() - ((v.x() - beamSpot->x0()) * px + (v.y() - beamSpot->y0()) * py) / pt * pz / pt;
0169               if (fabs(z) < m_maxZ) {
0170                 zProjections.push_back(z);
0171               }
0172             }
0173           }  //if compatible cluster
0174         }  // loop on module hits
0175       }  // if compatible module
0176     }  // loop on pixel modules
0177
0178   }  // loop on selected jets
0179   std::sort(zProjections.begin(), zProjections.end());
0180
0181   std::vector<float>::iterator itCenter = zProjections.begin();
0182   std::vector<float>::iterator itLeftSide = zProjections.begin();
0183   std::vector<float>::iterator itRightSide = zProjections.begin();
0184   std::vector<int> counts;
0185   float zCluster = m_clusterLength / 2.0;  //cm
0186   int max = 0;
0187   std::vector<float>::iterator left, right;
0188   for (; itCenter != zProjections.end(); itCenter++) {
0189     while (itLeftSide != zProjections.end() && (*itCenter - *itLeftSide) > zCluster)
0190       itLeftSide++;
0191     while (itRightSide != zProjections.end() && (*itRightSide - *itCenter) < zCluster)
0192       itRightSide++;
0193
0194     int n = itRightSide - itLeftSide;
0195     // std::cout << "algo :"<< *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl;
0196     counts.push_back(n);
0197     if (n > max) {
0198       max = n;
0199       left = itLeftSide;
0200     }
0201     if (n >= max) {
0202       max = n;
0203       right = itRightSide;
0204       //          std::cout << "algo :"<< i << " " << j << " " << *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl;
0205     }
0206   }
0207
0208   float res = 0;
0209   if (!zProjections.empty()) {
0210     res = *(left + (right - left) / 2);
0211     //     std::cout << "RES " << res << std::endl;
0212     Vertex::Error e;
0213     e(0, 0) = 0.0015 * 0.0015;
0214     e(1, 1) = 0.0015 * 0.0015;
0215     e(2, 2) = 1.5 * 1.5;
0216     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
0217     Vertex thePV(p, e, 1, 1, 0);
0218     auto pOut = std::make_unique<reco::VertexCollection>();
0219     pOut->push_back(thePV);
0220     iEvent.put(std::move(pOut));
0221   } else {
0222     //   std::cout << "DUMMY " << res << std::endl;
0223
0224     Vertex::Error e;
0225     e(0, 0) = 0.0015 * 0.0015;
0226     e(1, 1) = 0.0015 * 0.0015;
0227     e(2, 2) = 1.5 * 1.5;
0228     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
0229     Vertex thePV(p, e, 0, 0, 0);
0230     auto pOut = std::make_unique<reco::VertexCollection>();
0231     pOut->push_back(thePV);
0232     iEvent.put(std::move(pOut));
0233   }
0234 }
0235
0236 //define this as a plug-in
0237 DEFINE_FWK_MODULE(FastPrimaryVertexProducer);