FastPrimaryVertexProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    FastPrimaryVertexProducer
// Class:      FastPrimaryVertexProducer
//
//
// Original Author:  Andrea RIZZI
//         Created:  Thu Dec 22 14:51:44 CET 2011
//
//


// system include files
#include <memory>
#include <vector>
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"

#include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"

#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
#include "DataFormats/Math/interface/deltaPhi.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"

#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"

#include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
#include "CommonTools/Clustering1D/interface/Cluster1DMerger.h"
#include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
#define HaveMtv
#define HaveFsmw
#define HaveDivisive
#ifdef HaveMtv
#include "CommonTools/Clustering1D/interface/MtvClusterizer1D.h"
#endif
#ifdef HaveFsmw
#include "CommonTools/Clustering1D/interface/FsmwClusterizer1D.h"
#endif
#ifdef HaveDivisive
#include "CommonTools/Clustering1D/interface/DivisiveClusterizer1D.h"
#endif

using namespace std;

//
// class declaration
//

class FastPrimaryVertexProducer : public edm::global::EDProducer<> {
   public:
      explicit FastPrimaryVertexProducer(const edm::ParameterSet&);

   private:
      void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
      edm::EDGetTokenT<SiPixelClusterCollectionNew> m_clusters;
      edm::EDGetTokenT<edm::View<reco::Jet> > m_jets;
      edm::EDGetTokenT<reco::BeamSpot> m_beamSpot;
      std::string m_pixelCPE;
      double m_maxZ;
      double m_maxSizeX;
      double m_maxDeltaPhi;
      double m_clusterLength;

};

FastPrimaryVertexProducer::FastPrimaryVertexProducer(const edm::ParameterSet& iConfig)
{
  m_clusters          = consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusters"));
  m_jets              = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
  m_beamSpot          = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"));
  m_pixelCPE          = iConfig.getParameter<std::string>("pixelCPE");
  m_maxZ          = iConfig.getParameter<double>("maxZ");
  m_maxSizeX          = iConfig.getParameter<double>("maxSizeX");
  m_maxDeltaPhi       = iConfig.getParameter<double>("maxDeltaPhi");
  m_clusterLength     = iConfig.getParameter<double>("clusterLength");
  produces<reco::VertexCollection>();
}



void
FastPrimaryVertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const
{
   using namespace edm;
   using namespace reco;
   using namespace std;

   Handle<SiPixelClusterCollectionNew> cH;
   iEvent.getByToken(m_clusters,cH);
   const SiPixelClusterCollectionNew & pixelClusters = *cH.product();

   Handle<edm::View<reco::Jet> > jH;
   iEvent.getByToken(m_jets,jH);
   const edm::View<reco::Jet> & jets = *jH.product();

   CaloJetCollection selectedJets;
   for(edm::View<reco::Jet>::const_iterator it = jets.begin() ; it != jets.end() ; it++)
   {
    if(it->pt() > 40 && fabs(it->eta()) < 1.6)
    {
      const CaloJet * ca = dynamic_cast<const CaloJet *>( &(*it));
      if(ca ==nullptr) abort();
      selectedJets.push_back(*ca);
//    std::cout << "Jet eta,phi,pt: "<< it->eta() << "," << it->phi() << "," << it->pt()   << std::endl;
    }
   }

   edm::ESHandle<PixelClusterParameterEstimator> pe;
   const PixelClusterParameterEstimator * pp ;
   iSetup.get<TkPixelCPERecord>().get(m_pixelCPE , pe );
   pp = pe.product();

   edm::Handle<BeamSpot> beamSpot;
   iEvent.getByToken(m_beamSpot,beamSpot);

   edm::ESHandle<TrackerGeometry> tracker;
   iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
   const TrackerGeometry * trackerGeometry = tracker.product();


   float lengthBmodule=6.66;//cm
   std::vector<float> zProjections;
   for(CaloJetCollection::const_iterator jit = selectedJets.begin() ; jit != selectedJets.end() ; jit++)
   {
     float px=jit->px();
     float py=jit->py();
     float pz=jit->pz();
     float pt=jit->pt();

     float jetZOverRho = jit->momentum().Z()/jit->momentum().Rho();
     int minSizeY = fabs(2.*jetZOverRho)-1;
     int  maxSizeY = fabs(2.*jetZOverRho)+2;
     if( fabs(jit->eta()) > 1.6)
     {
        minSizeY = 1;
     }

     for(SiPixelClusterCollectionNew::const_iterator it = pixelClusters.begin() ; it != pixelClusters.end() ; it++) //Loop on pixel modules with clusters
     {
        DetId id = it->detId();
        const edmNew::DetSet<SiPixelCluster> & detset  = (*it);
        Point3DBase<float, GlobalTag> modulepos=trackerGeometry->idToDet(id)->position();
        float zmodule = modulepos.z() - ((modulepos.x()-beamSpot->x0())*px+(modulepos.y()-beamSpot->y0())*py)/pt * pz/pt;
        if ((fabs(deltaPhi(jit->momentum().Phi(),modulepos.phi()))< m_maxDeltaPhi*2)&&(fabs(zmodule)<(m_maxZ+lengthBmodule/2))){

        for(size_t j = 0 ; j < detset.size() ; j ++) // Loop on pixel clusters on this module
        {
      const SiPixelCluster & aCluster =  detset[j];
          if(aCluster.sizeX() < m_maxSizeX && aCluster.sizeY() >= minSizeY && aCluster.sizeY() <= maxSizeY) {
            Point3DBase<float, GlobalTag> v = trackerGeometry->idToDet(id)->surface().toGlobal(pp->localParametersV( aCluster,( *trackerGeometry->idToDetUnit(id)))[0].first) ;
            GlobalPoint v_bs(v.x()-beamSpot->x0(),v.y()-beamSpot->y0(),v.z());
            if(fabs(deltaPhi(jit->momentum().Phi(),v_bs.phi())) < m_maxDeltaPhi)
            {
              float z = v.z() - ((v.x()-beamSpot->x0())*px+(v.y()-beamSpot->y0())*py)/pt * pz/pt;
              if(fabs(z) < m_maxZ)
              {
            zProjections.push_back(z);
          }
        }
          } //if compatible cluster
    } // loop on module hits
      } // if compatible module
    } // loop on pixel modules

   } // loop on selected jets
  std::sort(zProjections.begin(),zProjections.end());

  std::vector<float>::iterator itCenter = zProjections.begin();
  std::vector<float>::iterator itLeftSide = zProjections.begin();
  std::vector<float>::iterator itRightSide = zProjections.begin();
  std::vector<int> counts;
  float zCluster = m_clusterLength/2.0; //cm
  int max=0;
  std::vector<float>::iterator left,right;
  for(;itCenter!=zProjections.end(); itCenter++)
  {

    while(itLeftSide != zProjections.end() && (*itCenter - *itLeftSide) > zCluster  ) itLeftSide++;
    while(itRightSide != zProjections.end() && (*itRightSide - *itCenter) < zCluster  ) itRightSide++;

    int n= itRightSide-itLeftSide;
   // std::cout << "algo :"<< *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl;
    counts.push_back(n);
    if(n > max) {
         max=n;
         left=itLeftSide;
    }
    if(n >= max) {
      max=n;
          right=itRightSide;
//          std::cout << "algo :"<< i << " " << j << " " << *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl;
    }
  }




  float res=0;
  if(!zProjections.empty())
  {
     res=*(left+(right-left)/2);
//     std::cout << "RES " << res << std::endl;
     Vertex::Error e;
     e(0, 0) = 0.0015 * 0.0015;
     e(1, 1) = 0.0015 * 0.0015;
     e(2, 2) = 1.5 * 1.5;
     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
     Vertex thePV(p, e, 1, 1, 0);
     auto pOut = std::make_unique<reco::VertexCollection>();
     pOut->push_back(thePV);
     iEvent.put(std::move(pOut));
   } else
   {
  //   std::cout << "DUMMY " << res << std::endl;

     Vertex::Error e;
     e(0, 0) = 0.0015 * 0.0015;
     e(1, 1) = 0.0015 * 0.0015;
     e(2, 2) = 1.5 * 1.5;
     Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
     Vertex thePV(p, e, 0, 0, 0);
     auto pOut = std::make_unique<reco::VertexCollection>();
     pOut->push_back(thePV);
     iEvent.put(std::move(pOut));

   }


}


//define this as a plug-in
DEFINE_FWK_MODULE(FastPrimaryVertexProducer);