Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    FastPrimaryVertexProducer
00004 // Class:      FastPrimaryVertexProducer
00005 // 
00013 //
00014 // Original Author:  Andrea RIZZI
00015 //         Created:  Thu Dec 22 14:51:44 CET 2011
00016 // $Id: FastPrimaryVertexProducer.cc,v 1.3 2012/02/02 09:04:04 arizzi Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 #include <vector>
00024 // user include files
00025 #include "FWCore/Framework/interface/Frameworkfwd.h"
00026 #include "FWCore/Framework/interface/EDProducer.h"
00027 #include "FWCore/Framework/interface/ESHandle.h"
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/MakerMacros.h"
00030 
00031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00032 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
00033 #include "DataFormats/JetReco/interface/CaloJet.h"
00034 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00035 
00036 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
00037 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00038 
00039 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00040 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
00041 #include "DataFormats/Math/interface/deltaPhi.h"
00042 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00043 #include "DataFormats/VertexReco/interface/Vertex.h"
00044 #include "DataFormats/TrackReco/interface/Track.h"
00045 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
00046 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00047 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00048 
00049 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00050 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00051 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
00052 
00053 #include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
00054 #include "CommonTools/Clustering1D/interface/Cluster1DMerger.h"
00055 #include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
00056 #define HaveMtv
00057 #define HaveFsmw
00058 #define HaveDivisive
00059 #ifdef HaveMtv
00060 #include "CommonTools/Clustering1D/interface/MtvClusterizer1D.h"
00061 #endif
00062 #ifdef HaveFsmw
00063 #include "CommonTools/Clustering1D/interface/FsmwClusterizer1D.h"
00064 #endif
00065 #ifdef HaveDivisive
00066 #include "CommonTools/Clustering1D/interface/DivisiveClusterizer1D.h"
00067 #endif
00068 
00069 using namespace std;
00070 
00071 //
00072 // class declaration
00073 //
00074 
00075 class FastPrimaryVertexProducer : public edm::EDProducer {
00076    public:
00077       explicit FastPrimaryVertexProducer(const edm::ParameterSet&);
00078 
00079    private:
00080       virtual void produce(edm::Event&, const edm::EventSetup&);
00081       edm::InputTag m_clusters;
00082       edm::InputTag m_jets;
00083       edm::InputTag m_beamSpot;
00084       std::string m_pixelCPE; 
00085       double m_maxZ; 
00086       double m_maxSizeX;
00087       double m_maxDeltaPhi;
00088       double m_clusterLength;   
00089 };
00090 
00091 FastPrimaryVertexProducer::FastPrimaryVertexProducer(const edm::ParameterSet& iConfig)
00092 {
00093   m_clusters          = iConfig.getParameter<edm::InputTag>("clusters");
00094   m_jets              = iConfig.getParameter<edm::InputTag>("jets");
00095   m_beamSpot          = iConfig.getParameter<edm::InputTag>("beamSpot");
00096   m_pixelCPE          = iConfig.getParameter<std::string>("pixelCPE");
00097   m_maxZ              = iConfig.getParameter<double>("maxZ");   
00098   m_maxSizeX          = iConfig.getParameter<double>("maxSizeX");       
00099   m_maxDeltaPhi       = iConfig.getParameter<double>("maxDeltaPhi");    
00100   m_clusterLength     = iConfig.getParameter<double>("clusterLength");  
00101   produces<reco::VertexCollection>();
00102 }
00103 
00104 
00105 
00106 void
00107 FastPrimaryVertexProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00108 {
00109    using namespace edm;
00110    using namespace reco;
00111    using namespace std;
00112 
00113    Handle<SiPixelClusterCollectionNew> cH;
00114    iEvent.getByLabel(m_clusters,cH);
00115    const SiPixelClusterCollectionNew & pixelClusters = *cH.product();
00116 
00117    Handle<edm::View<reco::Jet> > jH;
00118    iEvent.getByLabel(m_jets,jH);
00119    const edm::View<reco::Jet> & jets = *jH.product();
00120 
00121    CaloJetCollection selectedJets;
00122    for(edm::View<reco::Jet>::const_iterator it = jets.begin() ; it != jets.end() ; it++)
00123    {
00124     if(it->pt() > 40 && fabs(it->eta()) < 1.6)
00125     {
00126       const CaloJet * ca = dynamic_cast<const CaloJet *>( &(*it));
00127       if(ca ==0) abort();
00128       selectedJets.push_back(*ca);
00129 //    std::cout << "Jet eta,phi,pt: "<< it->eta() << "," << it->phi() << "," << it->pt()   << std::endl;
00130     } 
00131    }
00132   
00133    edm::ESHandle<PixelClusterParameterEstimator> pe; 
00134    const PixelClusterParameterEstimator * pp ;
00135    iSetup.get<TkPixelCPERecord>().get(m_pixelCPE , pe );  
00136    pp = pe.product();
00137 
00138    edm::Handle<BeamSpot> beamSpot;
00139    iEvent.getByLabel(m_beamSpot,beamSpot);
00140  
00141    edm::ESHandle<TrackerGeometry> tracker;
00142    iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
00143    const TrackerGeometry * trackerGeometry = tracker.product();
00144 
00145 
00146    float lengthBmodule=6.66;//cm 
00147    std::vector<float> zProjections;
00148    for(CaloJetCollection::const_iterator jit = selectedJets.begin() ; jit != selectedJets.end() ; jit++)
00149    {
00150      float px=jit->px();
00151      float py=jit->py();
00152      float pz=jit->pz();
00153      float pt=jit->pt();
00154 
00155      float jetZOverRho = jit->momentum().Z()/jit->momentum().Rho();
00156      int minSizeY = fabs(2.*jetZOverRho)-1; 
00157      int  maxSizeY = fabs(2.*jetZOverRho)+2;
00158      if( fabs(jit->eta()) > 1.6)
00159      {
00160         minSizeY = 1;
00161      }
00162 
00163      for(SiPixelClusterCollectionNew::const_iterator it = pixelClusters.begin() ; it != pixelClusters.end() ; it++) //Loop on pixel modules with clusters
00164      {
00165         DetId id = it->detId();
00166         const edmNew::DetSet<SiPixelCluster> & detset  = (*it);
00167         Point3DBase<float, GlobalTag> modulepos=trackerGeometry->idToDet(id)->position();
00168         float zmodule = modulepos.z() - ((modulepos.x()-beamSpot->x0())*px+(modulepos.y()-beamSpot->y0())*py)/pt * pz/pt;
00169         if ((fabs(deltaPhi(jit->momentum().Phi(),modulepos.phi()))< m_maxDeltaPhi*2)&&(fabs(zmodule)<(m_maxZ+lengthBmodule/2))){
00170        
00171         for(size_t j = 0 ; j < detset.size() ; j ++) // Loop on pixel clusters on this module
00172         {
00173           const SiPixelCluster & aCluster =  detset[j];
00174           if(aCluster.sizeX() < m_maxSizeX && aCluster.sizeY() >= minSizeY && aCluster.sizeY() <= maxSizeY) {
00175             Point3DBase<float, GlobalTag> v = trackerGeometry->idToDet(id)->surface().toGlobal(pp->localParametersV( aCluster,( *trackerGeometry->idToDetUnit(id)))[0].first) ;
00176             GlobalPoint v_bs(v.x()-beamSpot->x0(),v.y()-beamSpot->y0(),v.z());
00177             if(fabs(deltaPhi(jit->momentum().Phi(),v_bs.phi())) < m_maxDeltaPhi)
00178             {
00179               float z = v.z() - ((v.x()-beamSpot->x0())*px+(v.y()-beamSpot->y0())*py)/pt * pz/pt;   
00180               if(fabs(z) < m_maxZ)
00181               {
00182                 zProjections.push_back(z); 
00183               } 
00184             }
00185           } //if compatible cluster   
00186         } // loop on module hits
00187       } // if compatible module
00188     } // loop on pixel modules
00189       
00190    } // loop on selected jets
00191   std::sort(zProjections.begin(),zProjections.end());
00192    
00193   std::vector<float>::iterator itCenter = zProjections.begin();
00194   std::vector<float>::iterator itLeftSide = zProjections.begin();
00195   std::vector<float>::iterator itRightSide = zProjections.begin();
00196   std::vector<int> counts;
00197   float zCluster = m_clusterLength/2.0; //cm 
00198   int max=0;
00199   std::vector<float>::iterator left,right;
00200   for(;itCenter!=zProjections.end(); itCenter++)
00201   {
00202   
00203     while(itLeftSide != zProjections.end() && (*itCenter - *itLeftSide) > zCluster  ) itLeftSide++;
00204     while(itRightSide != zProjections.end() && (*itRightSide - *itCenter) < zCluster  ) itRightSide++;
00205    
00206     int n= itRightSide-itLeftSide;
00207    // std::cout << "algo :"<< *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl; 
00208     counts.push_back(n);
00209     if(n > max) {
00210          max=n; 
00211          left=itLeftSide;
00212     }
00213     if(n >= max) {
00214           max=n;
00215           right=itRightSide;
00216 //          std::cout << "algo :"<< i << " " << j << " " << *itCenter << " " << itCenter-zProjections.begin() << "  dists: " <<  (*itCenter - *itLeftSide) << " " << (*itRightSide - *itCenter) << " count: " <<  n << std::endl;
00217     }
00218   }
00219 
00220 
00221  
00222 
00223   float res=0; 
00224   if(zProjections.size() > 0) 
00225   {
00226      res=*(left+(right-left)/2);
00227 //     std::cout << "RES " << res << std::endl;
00228      Vertex::Error e; 
00229      e(0, 0) = 0.0015 * 0.0015;
00230      e(1, 1) = 0.0015 * 0.0015;
00231      e(2, 2) = 1.5 * 1.5;
00232      Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
00233      Vertex thePV(p, e, 1, 1, 0);
00234      std::auto_ptr<reco::VertexCollection> pOut(new reco::VertexCollection());
00235      pOut->push_back(thePV);
00236      iEvent.put(pOut);
00237    } else
00238    {
00239   //   std::cout << "DUMMY " << res << std::endl;
00240 
00241      Vertex::Error e;
00242      e(0, 0) = 0.0015 * 0.0015;
00243      e(1, 1) = 0.0015 * 0.0015;
00244      e(2, 2) = 1.5 * 1.5;
00245      Vertex::Point p(beamSpot->x(res), beamSpot->y(res), res);
00246      Vertex thePV(p, e, 0, 0, 0);
00247      std::auto_ptr<reco::VertexCollection> pOut(new reco::VertexCollection());
00248      pOut->push_back(thePV);
00249      iEvent.put(pOut);
00250 
00251    }
00252 
00253  
00254 }
00255 
00256 
00257 //define this as a plug-in
00258 DEFINE_FWK_MODULE(FastPrimaryVertexProducer);
/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoPixelVertexing/PixelVertexFinding/src/FastPrimaryVertexProducer.cc
