CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MCvsRecoVerticesAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: MCvsRecoVerticesAnalyzer
4 // Class: MCvsRecoVerticesAnalyzer
5 //
13 //
14 // Original Author: Andrea Venturi
15 // Created: Thu Dec 16 16:32:56 CEST 2010
16 //
17 //
18 
19 
20 // system include files
21 #include <memory>
22 #include <numeric>
23 
24 #include <vector>
25 
26 // user include files
29 
33 
35 
37 
40 
42 
46 
49 
50 #include "TH1F.h"
51 #include "TH2F.h"
52 #include "TProfile.h"
53 
54 //
55 // class decleration
56 //
57 
58 
60 public:
63 
64 private:
65  virtual void beginJob() ;
66  virtual void beginRun(const edm::Run&, const edm::EventSetup&);
67  virtual void endRun(const edm::Run&, const edm::EventSetup&);
68  virtual void analyze(const edm::Event&, const edm::EventSetup&);
69  virtual void endJob() ;
70 
71  // ----------member data ---------------------------
72 
73 
74 
75  const bool m_useweight;
82 
92 
96 
97 
98 };
99 
100 //
101 // constants, enums and typedefs
102 //
103 
104 //
105 // static data member definitions
106 //
107 
108 //
109 // constructors and destructor
110 //
112  : m_useweight( iConfig.getParameter< bool >( "useWeight" ) )
113  , m_useVisibleVertices( iConfig.getParameter< bool >( "useVisibleVertices" ) )
114  , m_histoParameters( iConfig.getUntrackedParameter< edm::ParameterSet >( "histoParameters", edm::ParameterSet() ) )
115  , m_doubleToken( consumes< double >( iConfig.getParameter< edm::InputTag >( "weightProduct" ) ) )
116  , m_vecPileupSummaryInfoToken( consumes< std::vector<PileupSummaryInfo> >( iConfig.getParameter< edm::InputTag >( "pileupSummaryCollection" ) ) )
117  , m_recoVertexCollectionToken( consumes< reco::VertexCollection >( iConfig.getParameter< edm::InputTag >( "pvCollection" ) ) )
118  , m_hepMCProductToken( consumes< edm::HepMCProduct >( iConfig.getParameter< edm::InputTag >( "mcTruthCollection" ) ) )
119 {
120  //now do what ever initialization is needed
121 
122  if(m_useVisibleVertices) edm::LogInfo("UseVisibleVertices") << "Only visible vertices will be used to compute Npileup";
123 
125 
126  m_hrecovsmcnvtx2d = tfserv->make<TH2F>("recovsmcnvtx2d","Number of reco vertices vs pileup interactions",60,-0.5,59.5,60,-0.5,59.5);
127  m_hrecovsmcnvtx2d->GetXaxis()->SetTitle("Pileup Interactions"); m_hrecovsmcnvtx2d->GetYaxis()->SetTitle("Reco Vertices");
128  m_hrecovsmcnvtxprof = tfserv->make<TProfile>("recovsmcnvtxprof","Mean number of reco vs pileup vertices",60,-0.5,59.5);
129  m_hrecovsmcnvtxprof->GetXaxis()->SetTitle("Pileup Interactions"); m_hrecovsmcnvtxprof->GetYaxis()->SetTitle("Reco Vertices");
130 
131  m_hrecovsmclumi2d = tfserv->make<TH2F>("recovsmclumi2d","Number of reco vertices vs ave pileup interactions",200,0.,50.,60,-0.5,59.5);
132  m_hrecovsmclumi2d->GetXaxis()->SetTitle("Average Pileup Interactions"); m_hrecovsmclumi2d->GetYaxis()->SetTitle("Reco Vertices");
133  m_hrecovsmclumiprof = tfserv->make<TProfile>("recovsmclumiprof","Mean number of reco vs ave pileup vertices",200,0.,50.);
134  m_hrecovsmclumiprof->GetXaxis()->SetTitle("Average Pileup Interactions"); m_hrecovsmclumiprof->GetYaxis()->SetTitle("Reco Vertices");
135 
136  if(m_useweight) {
137  m_hrecovsmcnvtxweightedprof = tfserv->make<TProfile>("recovsmcnvtxweightedprof","Mean number of reco vs pileup vertices (1-w) weight",60,-0.5,59.5);
138  m_hrecovsmcnvtxweightedprof->GetXaxis()->SetTitle("Pileup Interactions"); m_hrecovsmcnvtxweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");
139 
140  m_hrecovsmclumiweightedprof = tfserv->make<TProfile>("recovsmclumiweightedprof","Mean number of reco vs ave pileup vertices (1-w) weight",
141  200,0.,50.);
142  m_hrecovsmclumiweightedprof->GetXaxis()->SetTitle("Average Pileup Interactions");
143  m_hrecovsmclumiweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");
144  }
145 
146  m_hdeltazfirst = tfserv->make<TH1F>("deltazfirst","Reco-MC vertex z position (first vertex)",
147  m_histoParameters.getUntrackedParameter<unsigned int>("zBins",1000),
148  m_histoParameters.getUntrackedParameter<double>("zMin",-1.),
149  m_histoParameters.getUntrackedParameter<double>("zMax",1.));
150  m_hdeltazfirst->GetXaxis()->SetTitle("#Delta z (cm)"); m_hdeltazfirst->GetYaxis()->SetTitle("Events");
151 
152  m_hdeltazclose = tfserv->make<TH1F>("deltazclose","Reco-MC vertex z position (closest vertex)",
153  m_histoParameters.getUntrackedParameter<unsigned int>("zBins",1000),
154  m_histoParameters.getUntrackedParameter<double>("zMin",-1.),
155  m_histoParameters.getUntrackedParameter<double>("zMax",1.));
156  m_hdeltazclose->GetXaxis()->SetTitle("#Delta z (cm)"); m_hdeltazclose->GetYaxis()->SetTitle("Events");
157 
158  m_hclosestvtx = tfserv->make<TH1F>("closestvtx","Closest reco vtx ID",30,-0.5,29.5);
159  m_hclosestvtx->GetXaxis()->SetTitle("Vtx ID"); m_hclosestvtx->GetYaxis()->SetTitle("Events");
160 
161  m_hdeltazfirstvsnpu = tfserv->make<TH2F>("deltazfirstvsnpu","Reco-MC vertex z position (first vertex) vs Npileup",30,-0.5,29.5,
162  m_histoParameters.getUntrackedParameter<unsigned int>("zBins",1000),
163  m_histoParameters.getUntrackedParameter<double>("zMin",-1.),
164  m_histoParameters.getUntrackedParameter<double>("zMax",1.));
165  m_hdeltazfirstvsnpu->GetXaxis()->SetTitle("pileup Interactions"); m_hdeltazfirstvsnpu->GetYaxis()->SetTitle("#Delta z (cm)");
166 
167  m_hdeltazclosevsnpu = tfserv->make<TH2F>("deltazclosevsnpu","Reco-MC vertex z position (closest vertex) v Npileup",30,-0.5,29.5,
168  m_histoParameters.getUntrackedParameter<unsigned int>("zBins",1000),
169  m_histoParameters.getUntrackedParameter<double>("zMin",-1.),
170  m_histoParameters.getUntrackedParameter<double>("zMax",1.));
171  m_hdeltazclosevsnpu->GetXaxis()->SetTitle("Pileup Interactions"); m_hdeltazclosevsnpu->GetYaxis()->SetTitle("#Delta z (cm)");
172 
173  m_hclosestvtxvsnpu = tfserv->make<TH2F>("closestvtxvsnpu","Closest reco vtx ID vs Npileup",30,-0.5,29.5,30,-0.5,29.5);
174  m_hclosestvtxvsnpu->GetXaxis()->SetTitle("Pileup Interactions"); m_hclosestvtxvsnpu->GetYaxis()->SetTitle("Vtx ID");
175 
176 }
177 
178 
180 {
181 
182  // do anything here that needs to be done at desctruction time
183  // (e.g. close files, deallocate resources etc.)
184 
185 }
186 
187 
188 //
189 // member functions
190 //
191 
192 // ------------ method called to for each event ------------
193 void
195 {
196 
197  double weight = 1.;
198 
199  if(m_useweight) {
200  edm::Handle< double > weightprod;
201  iEvent.getByToken( m_doubleToken, weightprod );
202 
203  weight = *weightprod;
204 
205  }
206 
208  iEvent.getByToken( m_vecPileupSummaryInfoToken, pileupinfos );
209 
210  // look for the intime PileupSummaryInfo
211 
212  std::vector<PileupSummaryInfo>::const_iterator pileupinfo;
213 
214  for(pileupinfo = pileupinfos->begin(); pileupinfo != pileupinfos->end() ; ++pileupinfo) {
215 
216  if(pileupinfo->getBunchCrossing()==0) break;
217 
218  }
219 
220  //
221 
223  iEvent.getByToken( m_recoVertexCollectionToken, pvcoll );
224 
225 
226  //
227 
228  if(pileupinfo->getBunchCrossing()!=0) {
229 
230  edm::LogError("NoInTimePileUpInfo") << "Cannot find the in-time pileup info " << pileupinfo->getBunchCrossing();
231 
232  }
233  else {
234 
235  int npileup = pileupinfo->getPU_NumInteractions();
236 
237  if(m_useVisibleVertices) npileup = pileupinfo->getPU_zpositions().size();
238 
239  m_hrecovsmcnvtx2d->Fill(npileup,pvcoll->size(),weight);
240  m_hrecovsmcnvtxprof->Fill(npileup,pvcoll->size(),weight);
241 
242  m_hrecovsmclumi2d->Fill(pileupinfo->getTrueNumInteractions(),pvcoll->size(),weight);
243  m_hrecovsmclumiprof->Fill(pileupinfo->getTrueNumInteractions(),pvcoll->size(),weight);
244 
245  if(m_useweight) {
246  m_hrecovsmcnvtxweightedprof->Fill(npileup,pvcoll->size(),1.-weight);
247  m_hrecovsmclumiweightedprof->Fill(pileupinfo->getTrueNumInteractions(),pvcoll->size(),1.-weight);
248  }
249  //
250 
252  iEvent.getByToken( m_hepMCProductToken, EvtHandle );
253 
254  const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
255 
256  // compute the difference between the main interaction vertex z position and the first vertex of the collection
257 
258  if(pvcoll->size() !=0) {
259  if(!(*pvcoll)[0].isFake()) {
260  // get the first vertex
261  if(Evt->vertices_begin() != Evt->vertices_end()) {
262  m_hdeltazfirst->Fill((*pvcoll)[0].z()-(*Evt->vertices_begin())->point3d().z()/10.,weight);
263  m_hdeltazfirstvsnpu->Fill(npileup,(*pvcoll)[0].z()-(*Evt->vertices_begin())->point3d().z()/10.,weight);
264  }
265  }
266  }
267 
268  // compute the difference between the main interaction vertex z position and the closest reco vertex
269 
270  double minabsdist = -1.;
271  double mindist = -999.;
272  int closestvtx = -1;
273 
274  for(unsigned int ivtx = 0 ; ivtx < pvcoll->size() ; ++ivtx) {
275 
276  if(closestvtx < 0 || minabsdist > std::abs((*pvcoll)[ivtx].z()-(*Evt->vertices_begin())->point3d().z()/10.)) {
277  mindist = (*pvcoll)[ivtx].z()-(*Evt->vertices_begin())->point3d().z()/10.;
278  closestvtx = ivtx;
279  minabsdist = std::abs(mindist);
280  }
281 
282  }
283  if(closestvtx >= 0) {
284  m_hdeltazclose->Fill(mindist,weight);
285  m_hdeltazclosevsnpu->Fill(npileup,mindist,weight);
286  m_hclosestvtx->Fill(closestvtx,weight);
287  m_hclosestvtxvsnpu->Fill(npileup,closestvtx,weight);
288  }
289 
290  }
291 }
292 
293  void
295 {
296 }
297 
298 void
300 {
301 }
302 
303 
304 
305 // ------------ method called once each job just before starting event loop ------------
306 void
308 {
309 }
310 
311 // ------------ method called once each job just after ending the event loop ------------
312 void
314 {
315 }
316 
317 //define this as a plug-in
T getUntrackedParameter(std::string const &, T const &) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
float float float z
edm::EDGetTokenT< reco::VertexCollection > m_recoVertexCollectionToken
int iEvent
Definition: GenABIO.cc:243
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
edm::EDGetTokenT< edm::HepMCProduct > m_hepMCProductToken
MCvsRecoVerticesAnalyzer(const edm::ParameterSet &)
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > m_vecPileupSummaryInfoToken
edm::EDGetTokenT< double > m_doubleToken
const edm::ParameterSet m_histoParameters
int weight
Definition: histoStyle.py:50
virtual void endRun(const edm::Run &, const edm::EventSetup &)
Definition: Run.h:41