test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PrimaryVertexAnalyzer4PUSlimmed.cc
Go to the documentation of this file.
2 
7 
8 // reco track and vertex
11 
12 // TrackingParticle
15 
16 // associator
19 
20 // DQM
22 
23 #include <numeric>
24 
25 namespace {
26  template <typename T, size_t N>
27  std::array<T, N+1> makeLogBins(const double min, const double max) {
28  const double minLog10 = std::log10(min);
29  const double maxLog10 = std::log10(max);
30  const double width = (maxLog10-minLog10)/N;
31  std::array<T, N+1> ret;
32  for(size_t i=0; i<= N; ++i) {
33  ret[i] = std::pow(10, minLog10 + i*width);
34  }
35  return ret;
36  }
37 }
38 
39 //
40 // constructors and destructor
41 //
43  const edm::ParameterSet& iConfig)
44  : verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
45  use_only_charged_tracks_(iConfig.getUntrackedParameter<bool>(
46  "use_only_charged_tracks", true)),
47  do_generic_sim_plots_(iConfig.getUntrackedParameter<bool>("do_generic_sim_plots")),
48  root_folder_(
49  iConfig.getUntrackedParameter<std::string>("root_folder",
50  "Validation/Vertices")),
51  vecPileupSummaryInfoToken_(consumes<std::vector<PileupSummaryInfo> >(
52  edm::InputTag(std::string("addPileupInfo")))),
53  trackingParticleCollectionToken_(consumes<TrackingParticleCollection>(
54  iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleCollection"))),
55  trackingVertexCollectionToken_(consumes<TrackingVertexCollection>(
56  iConfig.getUntrackedParameter<edm::InputTag>("trackingVertexCollection"))),
57  simToRecoAssociationToken_(consumes<reco::SimToRecoCollection>(
58  iConfig.getUntrackedParameter<edm::InputTag>("trackAssociatorMap"))),
59  recoToSimAssociationToken_(consumes<reco::RecoToSimCollection>(
60  iConfig.getUntrackedParameter<edm::InputTag>("trackAssociatorMap"))),
61  vertexAssociatorToken_(consumes<reco::VertexToTrackingVertexAssociator>(
62  iConfig.getUntrackedParameter<edm::InputTag>("vertexAssociator"))) {
63  reco_vertex_collections_ = iConfig.getParameter<std::vector<edm::InputTag> >(
64  "vertexRecoCollections");
65  for (auto const& l : reco_vertex_collections_) {
69  }
70 }
71 
73 
74 //
75 // member functions
76 //
78  DQMStore::IBooker& i, edm::Run const& iRun, edm::EventSetup const& iSetup) {
79  // TODO(rovere) make this booking method shorter and smarter,
80  // factorizing similar histograms with different prefix in a single
81  // method call.
82  float log_bins[31] = {
83  0.0, 0.0002, 0.0004, 0.0006, 0.0008, 0.001, 0.002,
84  0.004, 0.006, 0.008, 0.01, 0.02,
85  0.04, 0.06, 0.08, 0.1, 0.2,
86  0.4, 0.6, 0.8, 1.0, 2.0,
87  4.0, 6.0, 8.0, 10.0, 20.0,
88  40.0, 60.0, 80.0, 100.0
89  };
90  float log_mergez_bins[18] = {
91  0.0, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1,
92  0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0
93  };
94 
95  auto const log_pt_bins = makeLogBins<float, 100>(0.1, 1e4);
96  float log_pt2_bins[16] = {
97  0.0, 0.1, 0.5,
98  1.0, 2.0, 5.0,
99  10.0, 20.0, 50.0,
100  100.0, 200.0, 500.0,
101  1000.0, 2000.0, 5000.0,10000.0
102  };
103  float log_ntrk_bins[25] = {
104  0., 2.0, 4.0, 6.0, 8.0, 10.,
105  12.0, 14.0, 16.0, 18.0, 22.0,
106  26.0, 30.0, 35.0, 40.0,
107  45.0, 50.0, 55.0, 60.0, 70.0,
108  80.0, 90.0, 100.0, 150.0, 200.0
109  };
110 
111  // TODO(rovere) Possibly change or add the main DQMStore booking
112  // interface to allow booking a TProfile with variable bin-width
113  // using an array of floats, as done for the TH1F case, not of
114  // doubles.
115  double log_pt2_bins_double[16] = {
116  0.0, 0.1, 0.5,
117  1.0, 2.0, 5.0,
118  10.0, 20.0, 50.0,
119  100.0, 200.0, 500.0,
120  1000.0, 2000.0, 5000.0,10000.0
121  };
122 
125  mes_["root_folder"]["GenVtx_vs_BX"] =
126  i.book2D("GenVtx_vs_BX", "GenVtx_vs_BX", 16, -12.5, 3.5, 200, 0., 200.);
127  // Generated Primary Vertex Plots
128  mes_["root_folder"]["GenPV_X"] =
129  i.book1D("GenPV_X", "GeneratedPV_X", 120, -0.6, 0.6);
130  mes_["root_folder"]["GenPV_Y"] =
131  i.book1D("GenPV_Y", "GeneratedPV_Y", 120, -0.6, 0.6);
132  mes_["root_folder"]["GenPV_Z"] =
133  i.book1D("GenPV_Z", "GeneratedPV_Z", 120, -60., 60.);
134  mes_["root_folder"]["GenPV_R"] =
135  i.book1D("GenPV_R", "GeneratedPV_R", 120, 0, 0.6);
136  mes_["root_folder"]["GenPV_Pt2"] =
137  i.book1D("GenPV_Pt2", "GeneratedPV_Sum-pt2", 15, &log_pt2_bins[0]);
138  mes_["root_folder"]["GenPV_NumTracks"] =
139  i.book1D("GenPV_NumTracks", "GeneratedPV_NumTracks", 24, &log_ntrk_bins[0]);
140  mes_["root_folder"]["GenPV_ClosestDistanceZ"] =
141  i.book1D("GenPV_ClosestDistanceZ", "GeneratedPV_ClosestDistanceZ", 30,
142  &log_bins[0]);
143 
144  // All Generated Vertices, used for efficiency plots
145  mes_["root_folder"]["GenAllV_NumVertices"] = i.book1D(
146  "GenAllV_NumVertices", "GeneratedAllV_NumVertices", 100, 0., 200.);
147  mes_["root_folder"]["GenAllV_X"] =
148  i.book1D("GenAllV_X", "GeneratedAllV_X", 120, -0.6, 0.6);
149  mes_["root_folder"]["GenAllV_Y"] =
150  i.book1D("GenAllV_Y", "GeneratedAllV_Y", 120, -0.6, 0.6);
151  mes_["root_folder"]["GenAllV_Z"] =
152  i.book1D("GenAllV_Z", "GeneratedAllV_Z", 120, -60, 60);
153  mes_["root_folder"]["GenAllV_R"] =
154  i.book1D("GenAllV_R", "GeneratedAllV_R", 120, 0, 0.6);
155  mes_["root_folder"]["GenAllV_Pt2"] =
156  i.book1D("GenAllV_Pt2", "GeneratedAllV_Sum-pt2", 15, &log_pt2_bins[0]);
157  mes_["root_folder"]["GenAllV_NumTracks"] =
158  i.book1D("GenAllV_NumTracks", "GeneratedAllV_NumTracks", 24, &log_ntrk_bins[0]);
159  mes_["root_folder"]["GenAllV_ClosestDistanceZ"] =
160  i.book1D("GenAllV_ClosestDistanceZ", "GeneratedAllV_ClosestDistanceZ", 30,
161  &log_bins[0]);
162  mes_["root_folder"]["GenAllV_PairDistanceZ"] =
163  i.book1D("GenAllV_PairDistanceZ", "GeneratedAllV_PairDistanceZ",
164  1000, 0, 20);
165  mes_["root_folder"]["SignalIsHighestPt2"] =
166  i.book1D("SignalIsHighestPt2", "SignalIsHighestPt2", 2, -0.5, 1.5);
167  }
168 
169  for (auto const& l : reco_vertex_collections_) {
170  std::string label = l.label();
171  std::string current_folder = root_folder_ + "/" + label;
172  i.setCurrentFolder(current_folder);
173 
174  auto book1d = [&](const char *name, int bins, double min, double max) {
175  mes_[label][name] = i.book1D(name, name, bins, min, max);
176  };
177  auto book1dlogx = [&](const char *name, int bins, const float *xbinedges) {
178  mes_[label][name] = i.book1D(name, name, bins, xbinedges);
179  };
180  auto book2d = [&](const char *name,
181  int xbins, double xmin, double xmax,
182  int ybins, double ymin, double ymax) {
183  mes_[label][name] = i.book2D(name, name, xbins,xmin,xmax, ybins,ymin,ymax);
184  };
185  auto book2dlogx = [&](const char *name,
186  int xbins, const float *xbinedges,
187  int ybins, double ymin, double ymax) {
188  auto me = i.book2D(name, name, xbins,xbinedges[0],xbinedges[xbins], ybins,ymin,ymax);
189  me->getTH2F()->GetXaxis()->Set(xbins, xbinedges);
190  mes_[label][name] = me;
191  };
192 
193 
194  mes_[label]["RecoVtx_vs_GenVtx"] = i.bookProfile(
195  "RecoVtx_vs_GenVtx", "RecoVtx_vs_GenVtx", 125, 0., 250., 250, 0., 250.);
196  mes_[label]["MatchedRecoVtx_vs_GenVtx"] =
197  i.bookProfile("MatchedRecoVtx_vs_GenVtx", "MatchedRecoVtx_vs_GenVtx",
198  125, 0., 250., 250, 0., 250.);
199  mes_[label]["KindOfSignalPV"] =
200  i.book1D("KindOfSignalPV", "KindOfSignalPV", 9, -0.5, 8.5);
201  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(1, "!Highest!Assoc2Any");
202  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(2, "Highest!Assoc2Any");
203  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(3, "!HighestAssoc2First");
204  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(4, "HighestAssoc2First");
205  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(5, "!HighestAssoc2!First");
206  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(6, "HighestAssoc2!First");
207  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(7, "!HighestAssoc2First");
208  mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(8, "HighestAssoc2First");
209  mes_[label]["MisTagRate"] =
210  i.book1D("MisTagRate", "MisTagRate", 2, -0.5, 1.5);
211  mes_[label]["MisTagRate_vs_PU"] =
212  i.bookProfile("MisTagRate_vs_PU", "MisTagRate_vs_PU", 125, 0., 250.,
213  2, 0., 1.);
214  mes_[label]["MisTagRate_vs_sum-pt2"] =
215  i.bookProfile("MisTagRate_vs_sum-pt2", "MisTagRate_vs_sum-pt2",
216  15, &log_pt2_bins_double[0], 2, 0., 1.);
217  mes_[label]["MisTagRate_vs_Z"] =
218  i.bookProfile("MisTagRate_vs_Z", "MisTagRate_vs_Z",
219  120, -60., 60., 2, 0., 1.);
220  mes_[label]["MisTagRate_vs_R"] =
221  i.bookProfile("MisTagRate_vs_R", "MisTagRate_vs_R",
222  120, 0., 0.6, 2, 0., 1.);
223  mes_[label]["MisTagRate_vs_NumTracks"] =
224  i.bookProfile("MisTagRate_vs_NumTracks", "MisTagRate_vs_NumTracks",
225  100, 0., 200, 2, 0., 1.);
226  mes_[label]["MisTagRateSignalIsHighest"] =
227  i.book1D("MisTagRateSignalIsHighest",
228  "MisTagRateSignalIsHighest", 2, -0.5, 1.5);
229  mes_[label]["MisTagRateSignalIsHighest_vs_PU"] =
230  i.bookProfile("MisTagRateSignalIsHighest_vs_PU",
231  "MisTagRateSignalIsHighest_vs_PU", 125, 0., 250.,
232  2, 0., 1.);
233  mes_[label]["MisTagRateSignalIsHighest_vs_sum-pt2"] =
234  i.bookProfile("MisTagRateSignalIsHighest_vs_sum-pt2",
235  "MisTagRateSignalIsHighest_vs_sum-pt2",
236  15, &log_pt2_bins_double[0], 2, 0., 1.);
237  mes_[label]["MisTagRateSignalIsHighest_vs_Z"] =
238  i.bookProfile("MisTagRateSignalIsHighest_vs_Z",
239  "MisTagRateSignalIsHighest_vs_Z",
240  120, -60., 60., 2, 0., 1.);
241  mes_[label]["MisTagRateSignalIsHighest_vs_R"] =
242  i.bookProfile("MisTagRateSignalIsHighest_vs_R",
243  "MisTagRateSignalIsHighest_vs_R",
244  120, 0., 0.6, 2, 0., 1.);
245  mes_[label]["MisTagRateSignalIsHighest_vs_NumTracks"] =
246  i.bookProfile("MisTagRateSignalIsHighest_vs_NumTracks",
247  "MisTagRateSignalIsHighest_vs_NumTracks",
248  100, 0., 200, 2, 0., 1.);
249  mes_[label]["MisTagRateSignalIsNotHighest"] =
250  i.book1D("MisTagRateSignalIsNotHighest",
251  "MisTagRateSignalIsNotHighest", 2, -0.5, 1.5);
252  mes_[label]["MisTagRateSignalIsNotHighest_vs_PU"] =
253  i.bookProfile("MisTagRateSignalIsNotHighest_vs_PU",
254  "MisTagRateSignalIsNotHighest_vs_PU", 125, 0., 250.,
255  2, 0., 1.);
256  mes_[label]["MisTagRateSignalIsNotHighest_vs_sum-pt2"] =
257  i.bookProfile("MisTagRateSignalIsNotHighest_vs_sum-pt2",
258  "MisTagRateSignalIsNotHighest_vs_sum-pt2",
259  15, &log_pt2_bins_double[0], 2, 0., 1.);
260  mes_[label]["MisTagRateSignalIsNotHighest_vs_Z"] =
261  i.bookProfile("MisTagRateSignalIsNotHighest_vs_Z",
262  "MisTagRateSignalIsNotHighest_vs_Z",
263  120, -60., 60., 2, 0., 1.);
264  mes_[label]["MisTagRateSignalIsNotHighest_vs_R"] =
265  i.bookProfile("MisTagRateSignalIsNotHighest_vs_R",
266  "MisTagRateSignalIsNotHighest_vs_R",
267  120, 0., 0.6, 2, 0., 1.);
268  mes_[label]["MisTagRateSignalIsNotHighest_vs_NumTracks"] =
269  i.bookProfile("MisTagRateSignalIsNotHighest_vs_NumTracks",
270  "MisTagRateSignalIsNotHighest_vs_NumTracks",
271  100, 0., 200, 2, 0., 1.);
272  mes_[label]["TruePVLocationIndex"] =
273  i.book1D("TruePVLocationIndex",
274  "TruePVLocationIndexInRecoVertexCollection", 12, -1.5, 10.5);
275  mes_[label]["TruePVLocationIndexCumulative"] =
276  i.book1D("TruePVLocationIndexCumulative",
277  "TruePVLocationIndexInRecoVertexCollectionCumulative",
278  3, -1.5, 1.5);
279  mes_[label]["TruePVLocationIndexSignalIsHighest"] =
280  i.book1D("TruePVLocationIndexSignalIsHighest",
281  "TruePVLocationIndexSignalIsHighestInRecoVertexCollection",
282  12, -1.5, 10.5);
283  mes_[label]["TruePVLocationIndexSignalIsNotHighest"] =
284  i.book1D("TruePVLocationIndexSignalIsNotHighest",
285  "TruePVLocationIndexSignalIsNotHighestInRecoVertexCollection",
286  12, -1.5, 10.5);
287  // All Generated Vertices. Used for Efficiency plots We kind of
288  // duplicate plots here in case we want to perform more detailed
289  // studies on a selection of generated vertices, not on all of them.
290  mes_[label]["GenAllAssoc2Reco_NumVertices"] =
291  i.book1D("GenAllAssoc2Reco_NumVertices",
292  "GeneratedAllAssoc2Reco_NumVertices", 100, 0., 200.);
293  mes_[label]["GenAllAssoc2Reco_X"] = i.book1D(
294  "GenAllAssoc2Reco_X", "GeneratedAllAssoc2Reco_X", 120, -0.6, 0.6);
295  mes_[label]["GenAllAssoc2Reco_Y"] = i.book1D(
296  "GenAllAssoc2Reco_Y", "GeneratedAllAssoc2Reco_Y", 120, -0.6, 0.6);
297  mes_[label]["GenAllAssoc2Reco_Z"] = i.book1D(
298  "GenAllAssoc2Reco_Z", "GeneratedAllAssoc2Reco_Z", 120, -60, 60);
299  mes_[label]["GenAllAssoc2Reco_R"] =
300  i.book1D("GenAllAssoc2Reco_R", "GeneratedAllAssoc2Reco_R", 120, 0, 0.6);
301  mes_[label]["GenAllAssoc2Reco_Pt2"] =
302  i.book1D("GenAllAssoc2Reco_Pt2", "GeneratedAllAssoc2Reco_Sum-pt2", 15,
303  &log_pt2_bins[0]);
304  mes_[label]["GenAllAssoc2Reco_NumTracks"] =
305  i.book1D("GenAllAssoc2Reco_NumTracks",
306  "GeneratedAllAssoc2Reco_NumTracks", 24, &log_ntrk_bins[0]);
307  mes_[label]["GenAllAssoc2Reco_ClosestDistanceZ"] =
308  i.book1D("GenAllAssoc2Reco_ClosestDistanceZ",
309  "GeneratedAllAssoc2Reco_ClosestDistanceZ", 30, &log_bins[0]);
310  book1d("GenPVAssoc2RecoPV_X", 120, -0.6, 0.6);
311  book1d("GenPVAssoc2RecoPV_Y", 120, -0.6, 0.6);
312  book1d("GenPVAssoc2RecoPV_Z", 120, -60, 60);
313 
314  // All Generated Vertices Matched to a Reconstructed vertex. Used
315  // for Efficiency plots
316  mes_[label]["GenAllAssoc2RecoMatched_NumVertices"] =
317  i.book1D("GenAllAssoc2RecoMatched_NumVertices",
318  "GeneratedAllAssoc2RecoMatched_NumVertices", 100, 0., 200.);
319  mes_[label]["GenAllAssoc2RecoMatched_X"] =
320  i.book1D("GenAllAssoc2RecoMatched_X", "GeneratedAllAssoc2RecoMatched_X",
321  120, -0.6, 0.6);
322  mes_[label]["GenAllAssoc2RecoMatched_Y"] =
323  i.book1D("GenAllAssoc2RecoMatched_Y", "GeneratedAllAssoc2RecoMatched_Y",
324  120, -0.6, 0.6);
325  mes_[label]["GenAllAssoc2RecoMatched_Z"] =
326  i.book1D("GenAllAssoc2RecoMatched_Z", "GeneratedAllAssoc2RecoMatched_Z",
327  120, -60, 60);
328  mes_[label]["GenAllAssoc2RecoMatched_R"] =
329  i.book1D("GenAllAssoc2RecoMatched_R", "GeneratedAllAssoc2RecoMatched_R",
330  120, 0, 0.6);
331  mes_[label]["GenAllAssoc2RecoMatched_Pt2"] =
332  i.book1D("GenAllAssoc2RecoMatched_Pt2",
333  "GeneratedAllAssoc2RecoMatched_Sum-pt2", 15, &log_pt2_bins[0]);
334  mes_[label]["GenAllAssoc2RecoMatched_NumTracks"] =
335  i.book1D("GenAllAssoc2RecoMatched_NumTracks",
336  "GeneratedAllAssoc2RecoMatched_NumTracks", 24, &log_ntrk_bins[0]);
337  mes_[label]["GenAllAssoc2RecoMatched_ClosestDistanceZ"] = i.book1D(
338  "GenAllAssoc2RecoMatched_ClosestDistanceZ",
339  "GeneratedAllAssoc2RecoMatched_ClosestDistanceZ", 30, &log_bins[0]);
340  book1d("GenPVAssoc2RecoPVMatched_X", 120, -0.6, 0.6);
341  book1d("GenPVAssoc2RecoPVMatched_Y", 120, -0.6, 0.6);
342  book1d("GenPVAssoc2RecoPVMatched_Z", 120, -60, 60);
343 
344  // All Generated Vertices Multi-Matched to a Reconstructed vertex. Used
345  // for Duplicate rate plots
346  mes_[label]["GenAllAssoc2RecoMultiMatched_NumVertices"] = i.book1D(
347  "GenAllAssoc2RecoMultiMatched_NumVertices",
348  "GeneratedAllAssoc2RecoMultiMatched_NumVertices", 100, 0., 200.);
349  mes_[label]["GenAllAssoc2RecoMultiMatched_X"] =
350  i.book1D("GenAllAssoc2RecoMultiMatched_X",
351  "GeneratedAllAssoc2RecoMultiMatched_X", 120, -0.6, 0.6);
352  mes_[label]["GenAllAssoc2RecoMultiMatched_Y"] =
353  i.book1D("GenAllAssoc2RecoMultiMatched_Y",
354  "GeneratedAllAssoc2RecoMultiMatched_Y", 120, -0.6, 0.6);
355  mes_[label]["GenAllAssoc2RecoMultiMatched_Z"] =
356  i.book1D("GenAllAssoc2RecoMultiMatched_Z",
357  "GeneratedAllAssoc2RecoMultiMatched_Z", 120, -60, 60);
358  mes_[label]["GenAllAssoc2RecoMultiMatched_R"] =
359  i.book1D("GenAllAssoc2RecoMultiMatched_R",
360  "GeneratedAllAssoc2RecoMultiMatched_R", 120, 0, 0.6);
361  mes_[label]["GenAllAssoc2RecoMultiMatched_Pt2"] =
362  i.book1D("GenAllAssoc2RecoMultiMatched_Pt2",
363  "GeneratedAllAssoc2RecoMultiMatched_Sum-pt2",
364  15, &log_pt2_bins[0]);
365  mes_[label]["GenAllAssoc2RecoMultiMatched_NumTracks"] =
366  i.book1D("GenAllAssoc2RecoMultiMatched_NumTracks",
367  "GeneratedAllAssoc2RecoMultiMatched_NumTracks", 24, &log_ntrk_bins[0]);
368  mes_[label]["GenAllAssoc2RecoMultiMatched_ClosestDistanceZ"] = i.book1D(
369  "GenAllAssoc2RecoMultiMatched_ClosestDistanceZ",
370  "GeneratedAllAssoc2RecoMultiMatched_ClosestDistanceZ",
371  30, &log_bins[0]);
372 
373  // All Reco Vertices. Used for {Fake,Duplicate}-Rate plots
374  mes_[label]["RecoAllAssoc2Gen_NumVertices"] =
375  i.book1D("RecoAllAssoc2Gen_NumVertices",
376  "ReconstructedAllAssoc2Gen_NumVertices", 100, 0., 200.);
377  mes_[label]["RecoAllAssoc2Gen_X"] = i.book1D(
378  "RecoAllAssoc2Gen_X", "ReconstructedAllAssoc2Gen_X", 120, -0.6, 0.6);
379  mes_[label]["RecoAllAssoc2Gen_Y"] = i.book1D(
380  "RecoAllAssoc2Gen_Y", "ReconstructedAllAssoc2Gen_Y", 120, -0.6, 0.6);
381  mes_[label]["RecoAllAssoc2Gen_Z"] = i.book1D(
382  "RecoAllAssoc2Gen_Z", "ReconstructedAllAssoc2Gen_Z", 120, -60, 60);
383  mes_[label]["RecoAllAssoc2Gen_R"] = i.book1D(
384  "RecoAllAssoc2Gen_R", "ReconstructedAllAssoc2Gen_R", 120, 0, 0.6);
385  mes_[label]["RecoAllAssoc2Gen_Pt2"] =
386  i.book1D("RecoAllAssoc2Gen_Pt2", "ReconstructedAllAssoc2Gen_Sum-pt2",
387  15, &log_pt2_bins[0]);
388  mes_[label]["RecoAllAssoc2Gen_Ndof"] =
389  i.book1D("RecoAllAssoc2Gen_Ndof",
390  "ReconstructedAllAssoc2Gen_Ndof", 120, 0., 240.);
391  mes_[label]["RecoAllAssoc2Gen_NumTracks"] =
392  i.book1D("RecoAllAssoc2Gen_NumTracks",
393  "ReconstructedAllAssoc2Gen_NumTracks", 24, &log_ntrk_bins[0]);
394  mes_[label]["RecoAllAssoc2Gen_PU"] =
395  i.book1D("RecoAllAssoc2Gen_PU",
396  "ReconstructedAllAssoc2Gen_PU", 125, 0., 250.);
397  mes_[label]["RecoAllAssoc2Gen_ClosestDistanceZ"] =
398  i.book1D("RecoAllAssoc2Gen_ClosestDistanceZ",
399  "ReconstructedAllAssoc2Gen_ClosestDistanceZ",
400  30, &log_bins[0]);
401  mes_[label]["RecoAllAssoc2GenProperties"] =
402  i.book1D("RecoAllAssoc2GenProperties",
403  "ReconstructedAllAssoc2Gen_Properties", 8, -0.5, 7.5);
404  mes_[label]["RecoAllAssoc2Gen_PairDistanceZ"] =
405  i.book1D("RecoAllAssoc2Gen_PairDistanceZ",
406  "RecoAllAssoc2Gen_PairDistanceZ", 1000, 0, 20);
407 
408  // All Reconstructed Vertices Matched to a Generated vertex. Used
409  // for Fake-Rate plots
410  mes_[label]["RecoAllAssoc2GenMatched_NumVertices"] =
411  i.book1D("RecoAllAssoc2GenMatched_NumVertices",
412  "ReconstructedAllAssoc2GenMatched_NumVertices", 100, 0., 200.);
413  mes_[label]["RecoAllAssoc2GenMatched_X"] =
414  i.book1D("RecoAllAssoc2GenMatched_X",
415  "ReconstructedAllAssoc2GenMatched_X", 120, -0.6, 0.6);
416  mes_[label]["RecoAllAssoc2GenMatched_Y"] =
417  i.book1D("RecoAllAssoc2GenMatched_Y",
418  "ReconstructedAllAssoc2GenMatched_Y", 120, -0.6, 0.6);
419  mes_[label]["RecoAllAssoc2GenMatched_Z"] =
420  i.book1D("RecoAllAssoc2GenMatched_Z",
421  "ReconstructedAllAssoc2GenMatched_Z", 120, -60, 60);
422  mes_[label]["RecoAllAssoc2GenMatched_R"] =
423  i.book1D("RecoAllAssoc2GenMatched_R",
424  "ReconstructedAllAssoc2GenMatched_R", 120, 0, 0.6);
425  mes_[label]["RecoAllAssoc2GenMatched_Pt2"] =
426  i.book1D("RecoAllAssoc2GenMatched_Pt2",
427  "ReconstructedAllAssoc2GenMatched_Sum-pt2",
428  15, &log_pt2_bins[0]);
429  mes_[label]["RecoAllAssoc2GenMatched_Ndof"] =
430  i.book1D("RecoAllAssoc2GenMatched_Ndof",
431  "ReconstructedAllAssoc2GenMatched_Ndof", 120, 0., 240.);
432  mes_[label]["RecoAllAssoc2GenMatched_NumTracks"] =
433  i.book1D("RecoAllAssoc2GenMatched_NumTracks",
434  "ReconstructedAllAssoc2GenMatched_NumTracks", 24, &log_ntrk_bins[0]);
435  mes_[label]["RecoAllAssoc2GenMatched_PU"] =
436  i.book1D("RecoAllAssoc2GenMatched_PU",
437  "ReconstructedAllAssoc2GenMatched_PU", 125, 0., 250.);
438  mes_[label]["RecoAllAssoc2GenMatched_ClosestDistanceZ"] = i.book1D(
439  "RecoAllAssoc2GenMatched_ClosestDistanceZ",
440  "ReconstructedAllAssoc2GenMatched_ClosestDistanceZ", 30, &log_bins[0]);
441 
442  // All Reconstructed Vertices Multi-Matched to a Generated vertex. Used
443  // for Merge-Rate plots
444  mes_[label]["RecoAllAssoc2GenMultiMatched_NumVertices"] = i.book1D(
445  "RecoAllAssoc2GenMultiMatched_NumVertices",
446  "ReconstructedAllAssoc2GenMultiMatched_NumVertices", 100, 0., 200.);
447  mes_[label]["RecoAllAssoc2GenMultiMatched_X"] =
448  i.book1D("RecoAllAssoc2GenMultiMatched_X",
449  "ReconstructedAllAssoc2GenMultiMatched_X", 120, -0.6, 0.6);
450  mes_[label]["RecoAllAssoc2GenMultiMatched_Y"] =
451  i.book1D("RecoAllAssoc2GenMultiMatched_Y",
452  "ReconstructedAllAssoc2GenMultiMatched_Y", 120, -0.6, 0.6);
453  mes_[label]["RecoAllAssoc2GenMultiMatched_Z"] =
454  i.book1D("RecoAllAssoc2GenMultiMatched_Z",
455  "ReconstructedAllAssoc2GenMultiMatched_Z", 120, -60, 60);
456  mes_[label]["RecoAllAssoc2GenMultiMatched_R"] =
457  i.book1D("RecoAllAssoc2GenMultiMatched_R",
458  "ReconstructedAllAssoc2GenMultiMatched_R", 120, 0, 0.6);
459  mes_[label]["RecoAllAssoc2GenMultiMatched_Pt2"] = i.book1D(
460  "RecoAllAssoc2GenMultiMatched_Pt2",
461  "ReconstructedAllAssoc2GenMultiMatched_Sum-pt2", 15, &log_pt2_bins[0]);
462  mes_[label]["RecoAllAssoc2GenMultiMatched_NumTracks"] = i.book1D(
463  "RecoAllAssoc2GenMultiMatched_NumTracks",
464  "ReconstructedAllAssoc2GenMultiMatched_NumTracks", 24, &log_ntrk_bins[0]);
465  mes_[label]["RecoAllAssoc2GenMultiMatched_PU"] =
466  i.book1D("RecoAllAssoc2GenMultiMatched_PU",
467  "ReconstructedAllAssoc2GenMultiMatched_PU", 125, 0., 250.);
468  mes_[label]["RecoAllAssoc2GenMultiMatched_ClosestDistanceZ"] =
469  i.book1D("RecoAllAssoc2GenMultiMatched_ClosestDistanceZ",
470  "ReconstructedAllAssoc2GenMultiMatched_ClosestDistanceZ",
471  17, &log_mergez_bins[0]);
472 
473  // All Reconstructed Vertices Matched to a Multi-Matched Gen
474  // Vertex. Used for Duplicate rate plots done w.r.t. Reco
475  // Quantities. We basically want to ask how many times a RecoVTX
476  // has been reconstructed and associated to a SimulatedVTX that
477  // has been linked to at least another RecoVTX. In this sense this
478  // RecoVTX is a duplicate of the same, real GenVTX.
479  mes_[label]["RecoAllAssoc2MultiMatchedGen_NumVertices"] = i.book1D(
480  "RecoAllAssoc2MultiMatchedGen_NumVertices",
481  "RecoAllAssoc2MultiMatchedGen_NumVertices", 100, 0., 200.);
482  mes_[label]["RecoAllAssoc2MultiMatchedGen_X"] =
483  i.book1D("RecoAllAssoc2MultiMatchedGen_X",
484  "RecoAllAssoc2MultiMatchedGen_X", 120, -0.6, 0.6);
485  mes_[label]["RecoAllAssoc2MultiMatchedGen_Y"] =
486  i.book1D("RecoAllAssoc2MultiMatchedGen_Y",
487  "RecoAllAssoc2MultiMatchedGen_Y", 120, -0.6, 0.6);
488  mes_[label]["RecoAllAssoc2MultiMatchedGen_Z"] =
489  i.book1D("RecoAllAssoc2MultiMatchedGen_Z",
490  "RecoAllAssoc2MultiMatchedGen_Z", 120, -60, 60);
491  mes_[label]["RecoAllAssoc2MultiMatchedGen_R"] =
492  i.book1D("RecoAllAssoc2MultiMatchedGen_R",
493  "RecoAllAssoc2MultiMatchedGen_R", 120, 0, 0.6);
494  mes_[label]["RecoAllAssoc2MultiMatchedGen_Pt2"] =
495  i.book1D("RecoAllAssoc2MultiMatchedGen_Pt2",
496  "RecoAllAssoc2MultiMatchedGen_Sum-pt2", 15, &log_pt2_bins[0]);
497  mes_[label]["RecoAllAssoc2MultiMatchedGen_NumTracks"] =
498  i.book1D("RecoAllAssoc2MultiMatchedGen_NumTracks",
499  "RecoAllAssoc2MultiMatchedGen_NumTracks", 24, &log_ntrk_bins[0]);
500  mes_[label]["RecoAllAssoc2MultiMatchedGen_PU"] =
501  i.book1D("RecoAllAssoc2MultiMatchedGen_PU",
502  "RecoAllAssoc2MultiMatchedGen_PU", 125, 0., 250.);
503  mes_[label]["RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ"] = i.book1D(
504  "RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ",
505  "RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ", 30, &log_bins[0]);
506  mes_[label]["RecoAllAssoc2GenSimForMerge_ClosestDistanceZ"] = i.book1D(
507  "RecoAllAssoc2GenSimForMerge_ClosestDistanceZ",
508  "RecoAllAssoc2GenSimForMerge_ClosestDistanceZ",
509  17, &log_mergez_bins[0]);
510 
511 
512  // Resolution and pull histograms
513  const double resolx = 0.1;
514  const double resoly = 0.1;
515  const double resolz = 0.1;
516  const double resolpt2 = 10;
517 
518  // Non-merged vertices
519  book1d("RecoAllAssoc2GenMatched_ResolX", 100,-resolx,resolx);
520  book1d("RecoAllAssoc2GenMatched_ResolY", 100,-resoly,resoly);
521  book1d("RecoAllAssoc2GenMatched_ResolZ", 100,-resolz,resolz);
522  book1d("RecoAllAssoc2GenMatched_ResolPt2", 100,-resolpt2,resolpt2);
523 
524  book2d("RecoAllAssoc2GenMatched_ResolX_vs_PU", 125,0.,250., 100,-resolx,resolx);
525  book2d("RecoAllAssoc2GenMatched_ResolY_vs_PU", 125,0.,250., 100,-resoly,resoly);
526  book2d("RecoAllAssoc2GenMatched_ResolZ_vs_PU", 125,0.,250., 100,-resolz,resolz);
527  book2d("RecoAllAssoc2GenMatched_ResolPt2_vs_PU", 125,0.,250., 100,-resolpt2,resolpt2);
528 
529  book2dlogx("RecoAllAssoc2GenMatched_ResolX_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolx,resolx);
530  book2dlogx("RecoAllAssoc2GenMatched_ResolY_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resoly,resoly);
531  book2dlogx("RecoAllAssoc2GenMatched_ResolZ_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolz,resolz);
532  book2dlogx("RecoAllAssoc2GenMatched_ResolPt2_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolpt2,resolpt2);
533 
534  book2d("RecoAllAssoc2GenMatched_ResolX_vs_Z", 120,-60,60, 100,-resolx,resolx);
535  book2d("RecoAllAssoc2GenMatched_ResolY_vs_Z", 120,-60,60, 100,-resoly,resoly);
536  book2d("RecoAllAssoc2GenMatched_ResolZ_vs_Z", 120,-60,60, 100,-resolz,resolz);
537  book2d("RecoAllAssoc2GenMatched_ResolPt2_vs_Z", 120,-60,60, 100,-resolpt2,resolpt2);
538 
539  book2dlogx("RecoAllAssoc2GenMatched_ResolX_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolx,resolx);
540  book2dlogx("RecoAllAssoc2GenMatched_ResolY_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resoly,resoly);
541  book2dlogx("RecoAllAssoc2GenMatched_ResolZ_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolz,resolz);
542  book2dlogx("RecoAllAssoc2GenMatched_ResolPt2_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolpt2,resolpt2);
543 
544  book1d("RecoAllAssoc2GenMatched_PullX", 250,-25,25);
545  book1d("RecoAllAssoc2GenMatched_PullY", 250,-25,25);
546  book1d("RecoAllAssoc2GenMatched_PullZ", 250,-25,25);
547 
548  // Merged vertices
549  book1d("RecoAllAssoc2GenMatchedMerged_ResolX", 100,-resolx,resolx);
550  book1d("RecoAllAssoc2GenMatchedMerged_ResolY", 100,-resoly,resoly);
551  book1d("RecoAllAssoc2GenMatchedMerged_ResolZ", 100,-resolz,resolz);
552  book1d("RecoAllAssoc2GenMatchedMerged_ResolPt2", 100,-resolpt2,resolpt2);
553 
554  book2d("RecoAllAssoc2GenMatchedMerged_ResolX_vs_PU", 125,0.,250., 100,-resolx,resolx);
555  book2d("RecoAllAssoc2GenMatchedMerged_ResolY_vs_PU", 125,0.,250., 100,-resoly,resoly);
556  book2d("RecoAllAssoc2GenMatchedMerged_ResolZ_vs_PU", 125,0.,250., 100,-resolz,resolz);
557  book2d("RecoAllAssoc2GenMatchedMerged_ResolPt2_vs_PU", 125,0.,250., 100,-resolpt2,resolpt2);
558 
559  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolX_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolx,resolx);
560  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolY_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resoly,resoly);
561  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolZ_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolz,resolz);
562  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolPt2_vs_NumTracks", 24,&log_ntrk_bins[0], 100,-resolpt2,resolpt2);
563 
564  book2d("RecoAllAssoc2GenMatchedMerged_ResolX_vs_Z", 120,-60,60, 100,-resolx,resolx);
565  book2d("RecoAllAssoc2GenMatchedMerged_ResolY_vs_Z", 120,-60,60, 100,-resoly,resoly);
566  book2d("RecoAllAssoc2GenMatchedMerged_ResolZ_vs_Z", 120,-60,60, 100,-resolz,resolz);
567  book2d("RecoAllAssoc2GenMatchedMerged_ResolPt2_vs_Z", 120,-60,60, 100,-resolpt2,resolpt2);
568 
569  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolX_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolx,resolx);
570  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolY_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resoly,resoly);
571  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolZ_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolz,resolz);
572  book2dlogx("RecoAllAssoc2GenMatchedMerged_ResolPt2_vs_Pt", log_pt_bins.size()-1,&log_pt_bins[0], 100,-resolpt2,resolpt2);
573 
574  book1d("RecoAllAssoc2GenMatchedMerged_PullX", 250,-25,25);
575  book1d("RecoAllAssoc2GenMatchedMerged_PullY", 250,-25,25);
576  book1d("RecoAllAssoc2GenMatchedMerged_PullZ", 250,-25,25);
577 
578 
579  // Purity histograms
580  // Reco PV (vtx0) matched to hard-scatter gen vertex
581  book1d("RecoPVAssoc2GenPVMatched_Purity", 50, 0, 1);
582  book1d("RecoPVAssoc2GenPVMatched_Missing", 50, 0, 1);
583  book2d("RecoPVAssoc2GenPVMatched_Purity_vs_Index", 100,0,100, 50,0,1);
584 
585  // RECO PV (vtx0) not matched to hard-scatter gen vertex
586  book1d("RecoPVAssoc2GenPVNotMatched_Purity", 50, 0, 1);
587  book1d("RecoPVAssoc2GenPVNotMatched_Missing", 50, 0, 1);
588  book2d("RecoPVAssoc2GenPVNotMatched_Purity_vs_Index", 100,0,100, 50,0,1);
589 
590  // Purity vs. fake rate
591  book1d("RecoAllAssoc2Gen_Purity", 50, 0, 1); // denominator
592  book1d("RecoAllAssoc2GenMatched_Purity", 50, 0, 1); // 1-numerator
593 
594  // Vertex sum(pt2)
595  // The first two are orthogonal (i.e. their sum includes all reco vertices)
596  book1dlogx("RecoAssoc2GenPVMatched_Pt2", 15, &log_pt2_bins[0]);
597  book1dlogx("RecoAssoc2GenPVNotMatched_Pt2", 15, &log_pt2_bins[0]);
598 
599  book1dlogx("RecoAssoc2GenPVMatchedNotHighest_Pt2", 15, &log_pt2_bins[0]);
600  book1dlogx("RecoAssoc2GenPVNotMatched_GenPVTracksRemoved_Pt2", 15, &log_pt2_bins[0]);
601 
602  // Shared tracks
603  book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionReco", 50, 0, 1);
604  book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionReco", 50, 0, 1);
605  book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionRecoMatched", 50, 0, 1);
606  book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionRecoMatched", 50, 0, 1);
607 
608  book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionSim", 50, 0, 1);
609  book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionSim", 50, 0, 1);
610  book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionSimMatched", 50, 0, 1);
611  book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionSimMatched", 50, 0, 1);
612  }
613 }
614 
616  const simPrimaryVertex& v) {
617  if (v.eventId.event() == 0) {
618  mes_["root_folder"]["GenPV_X"]->Fill(v.x);
619  mes_["root_folder"]["GenPV_Y"]->Fill(v.y);
620  mes_["root_folder"]["GenPV_Z"]->Fill(v.z);
621  mes_["root_folder"]["GenPV_R"]->Fill(v.r);
622  mes_["root_folder"]["GenPV_Pt2"]->Fill(v.ptsq);
623  mes_["root_folder"]["GenPV_NumTracks"]->Fill(v.nGenTrk);
624  if (v.closest_vertex_distance_z > 0.)
625  mes_["root_folder"]["GenPV_ClosestDistanceZ"]
626  ->Fill(v.closest_vertex_distance_z);
627  }
628  mes_["root_folder"]["GenAllV_X"]->Fill(v.x);
629  mes_["root_folder"]["GenAllV_Y"]->Fill(v.y);
630  mes_["root_folder"]["GenAllV_Z"]->Fill(v.z);
631  mes_["root_folder"]["GenAllV_R"]->Fill(v.r);
632  mes_["root_folder"]["GenAllV_Pt2"]->Fill(v.ptsq);
633  mes_["root_folder"]["GenAllV_NumTracks"]->Fill(v.nGenTrk);
634  if (v.closest_vertex_distance_z > 0.)
635  mes_["root_folder"]["GenAllV_ClosestDistanceZ"]
636  ->Fill(v.closest_vertex_distance_z);
637 }
638 
640  const std::string& label,
642  mes_[label]["GenAllAssoc2Reco_X"]->Fill(v.x);
643  mes_[label]["GenAllAssoc2Reco_Y"]->Fill(v.y);
644  mes_[label]["GenAllAssoc2Reco_Z"]->Fill(v.z);
645  mes_[label]["GenAllAssoc2Reco_R"]->Fill(v.r);
646  mes_[label]["GenAllAssoc2Reco_Pt2"]->Fill(v.ptsq);
647  mes_[label]["GenAllAssoc2Reco_NumTracks"]->Fill(v.nGenTrk);
648  if (v.closest_vertex_distance_z > 0.)
649  mes_[label]["GenAllAssoc2Reco_ClosestDistanceZ"]
650  ->Fill(v.closest_vertex_distance_z);
651  if (v.rec_vertices.size()) {
652  mes_[label]["GenAllAssoc2RecoMatched_X"]->Fill(v.x);
653  mes_[label]["GenAllAssoc2RecoMatched_Y"]->Fill(v.y);
654  mes_[label]["GenAllAssoc2RecoMatched_Z"]->Fill(v.z);
655  mes_[label]["GenAllAssoc2RecoMatched_R"]->Fill(v.r);
656  mes_[label]["GenAllAssoc2RecoMatched_Pt2"]->Fill(v.ptsq);
657  mes_[label]["GenAllAssoc2RecoMatched_NumTracks"]->Fill(v.nGenTrk);
658  if (v.closest_vertex_distance_z > 0.)
659  mes_[label]["GenAllAssoc2RecoMatched_ClosestDistanceZ"]
660  ->Fill(v.closest_vertex_distance_z);
661  }
662  if (v.rec_vertices.size() > 1) {
663  mes_[label]["GenAllAssoc2RecoMultiMatched_X"]->Fill(v.x);
664  mes_[label]["GenAllAssoc2RecoMultiMatched_Y"]->Fill(v.y);
665  mes_[label]["GenAllAssoc2RecoMultiMatched_Z"]->Fill(v.z);
666  mes_[label]["GenAllAssoc2RecoMultiMatched_R"]->Fill(v.r);
667  mes_[label]["GenAllAssoc2RecoMultiMatched_Pt2"]->Fill(v.ptsq);
668  mes_[label]["GenAllAssoc2RecoMultiMatched_NumTracks"]->Fill(v.nGenTrk);
669  if (v.closest_vertex_distance_z > 0.)
670  mes_[label]["GenAllAssoc2RecoMultiMatched_ClosestDistanceZ"]
671  ->Fill(v.closest_vertex_distance_z);
672  }
673 }
674 
676  const std::string& label,
678  bool genPVMatchedToRecoPV) {
679  mes_[label]["GenPVAssoc2RecoPV_X"]->Fill(v.x);
680  mes_[label]["GenPVAssoc2RecoPV_Y"]->Fill(v.y);
681  mes_[label]["GenPVAssoc2RecoPV_Z"]->Fill(v.z);
682  if (genPVMatchedToRecoPV) {
683  mes_[label]["GenPVAssoc2RecoPVMatched_X"]->Fill(v.x);
684  mes_[label]["GenPVAssoc2RecoPVMatched_Y"]->Fill(v.y);
685  mes_[label]["GenPVAssoc2RecoPVMatched_Z"]->Fill(v.z);
686  }
687 }
688 
690  const std::string& label,
691  int num_pileup_vertices,
693  mes_[label]["RecoAllAssoc2Gen_X"]->Fill(v.x);
694  mes_[label]["RecoAllAssoc2Gen_Y"]->Fill(v.y);
695  mes_[label]["RecoAllAssoc2Gen_Z"]->Fill(v.z);
696  mes_[label]["RecoAllAssoc2Gen_R"]->Fill(v.r);
697  mes_[label]["RecoAllAssoc2Gen_Pt2"]->Fill(v.ptsq);
698  mes_[label]["RecoAllAssoc2Gen_Ndof"]->Fill(v.recVtx->ndof());
699  mes_[label]["RecoAllAssoc2Gen_NumTracks"]->Fill(v.nRecoTrk);
700  mes_[label]["RecoAllAssoc2Gen_PU"]->Fill(num_pileup_vertices);
701  mes_[label]["RecoAllAssoc2Gen_Purity"]->Fill(v.purity);
702  if (v.closest_vertex_distance_z > 0.)
703  mes_[label]["RecoAllAssoc2Gen_ClosestDistanceZ"]
704  ->Fill(v.closest_vertex_distance_z);
705  if (v.sim_vertices.size()) {
707  mes_[label]["RecoAllAssoc2GenMatched_X"]->Fill(v.x);
708  mes_[label]["RecoAllAssoc2GenMatched_Y"]->Fill(v.y);
709  mes_[label]["RecoAllAssoc2GenMatched_Z"]->Fill(v.z);
710  mes_[label]["RecoAllAssoc2GenMatched_R"]->Fill(v.r);
711  mes_[label]["RecoAllAssoc2GenMatched_Pt2"]->Fill(v.ptsq);
712  mes_[label]["RecoAllAssoc2GenMatched_Ndof"]->Fill(v.recVtx->ndof());
713  mes_[label]["RecoAllAssoc2GenMatched_NumTracks"]->Fill(v.nRecoTrk);
714  mes_[label]["RecoAllAssoc2GenMatched_PU"]->Fill(num_pileup_vertices);
715  mes_[label]["RecoAllAssoc2GenMatched_Purity"]->Fill(v.purity);
716  if (v.closest_vertex_distance_z > 0.)
717  mes_[label]["RecoAllAssoc2GenMatched_ClosestDistanceZ"]
718  ->Fill(v.closest_vertex_distance_z);
719 
720  // Fill resolution and pull plots here (as in MultiTrackValidator)
721  fillResolutionAndPullHistograms(label, num_pileup_vertices, v);
722 
723  // Now keep track of all RecoVTX associated to a SimVTX that
724  // itself is associated to more than one RecoVTX, for
725  // duplicate-rate plots on reco quantities.
726  if (v.sim_vertices_internal[0]->rec_vertices.size() > 1) {
728  mes_[label]["RecoAllAssoc2MultiMatchedGen_X"]->Fill(v.x);
729  mes_[label]["RecoAllAssoc2MultiMatchedGen_Y"]->Fill(v.y);
730  mes_[label]["RecoAllAssoc2MultiMatchedGen_Z"]->Fill(v.z);
731  mes_[label]["RecoAllAssoc2MultiMatchedGen_R"]->Fill(v.r);
732  mes_[label]["RecoAllAssoc2MultiMatchedGen_Pt2"]->Fill(v.ptsq);
733  mes_[label]["RecoAllAssoc2MultiMatchedGen_NumTracks"]->Fill(v.nRecoTrk);
734  mes_[label]["RecoAllAssoc2MultiMatchedGen_PU"]->Fill(num_pileup_vertices);
735  if (v.closest_vertex_distance_z > 0.)
736  mes_[label]["RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ"]
737  ->Fill(v.closest_vertex_distance_z);
738  }
739  // This is meant to be used as "denominator" for the merge-rate
740  // plots produced starting from reco quantities. We enter here
741  // only if the reco vertex has been associated, since we need info
742  // from the SimVTX associated to it. In this regard, the final
743  // merge-rate plot coming from reco is not to be intended as a
744  // pure efficiency-like plot, since the normalization is biased.
745  if (v.sim_vertices_internal[0]->closest_vertex_distance_z > 0.)
746  mes_[label]["RecoAllAssoc2GenSimForMerge_ClosestDistanceZ"]
747  ->Fill(v.sim_vertices_internal[0]->closest_vertex_distance_z);
748  }
749  // this plots are meant to be used to compute the merge rate
750  if (v.sim_vertices.size() > 1) {
752  mes_[label]["RecoAllAssoc2GenMultiMatched_X"]->Fill(v.x);
753  mes_[label]["RecoAllAssoc2GenMultiMatched_Y"]->Fill(v.y);
754  mes_[label]["RecoAllAssoc2GenMultiMatched_Z"]->Fill(v.z);
755  mes_[label]["RecoAllAssoc2GenMultiMatched_R"]->Fill(v.r);
756  mes_[label]["RecoAllAssoc2GenMultiMatched_Pt2"]->Fill(v.ptsq);
757  mes_[label]["RecoAllAssoc2GenMultiMatched_NumTracks"]->Fill(v.nRecoTrk);
758  mes_[label]["RecoAllAssoc2GenMultiMatched_PU"]->Fill(num_pileup_vertices);
759  if (v.sim_vertices_internal[0]->closest_vertex_distance_z > 0.)
760  mes_[label]["RecoAllAssoc2GenMultiMatched_ClosestDistanceZ"]
761  ->Fill(v.sim_vertices_internal[0]->closest_vertex_distance_z);
762  }
763  mes_[label]["RecoAllAssoc2GenProperties"]->Fill(v.kind_of_vertex);
764 
765 
766  std::string prefix;
767  if(v.sim_vertices.size() == 1) {
768  prefix = "RecoAllAssoc2GenSingleMatched_SharedTrackFraction";
769  }
770  else if(v.sim_vertices.size() > 1) {
771  prefix = "RecoAllAssoc2GenMultiMatched_SharedTrackFraction";
772  }
773 
774  for(size_t i=0; i<v.sim_vertices.size(); ++i) {
775  const double sharedTracks = v.sim_vertices_num_shared_tracks[i];
776  const simPrimaryVertex *simV = v.sim_vertices_internal[i];
777  mes_[label][prefix+"Reco"]->Fill(sharedTracks/v.nRecoTrk);
778  mes_[label][prefix+"RecoMatched"]->Fill(sharedTracks/v.num_matched_sim_tracks);
779  mes_[label][prefix+"Sim"]->Fill(sharedTracks/simV->nGenTrk);
780  mes_[label][prefix+"SimMatched"]->Fill(sharedTracks/simV->num_matched_reco_tracks);
781  }
782 }
783 
785  const std::string& label,
786  int num_pileup_vertices,
788 
789  std::string prefix = "RecoAllAssoc2GenMatched";
790  if(v.sim_vertices_internal.size() > 1) {
791  prefix += "Merged";
792  }
793 
794  // Use the best match as defined by the vertex truth associator
795  // reco-tracks as the best match
797  const double xres = v.x - bestMatch.x;
798  const double yres = v.y - bestMatch.y;
799  const double zres = v.z - bestMatch.z;
800  const double pt2res = v.ptsq - bestMatch.ptsq;
801 
802  const double xresol = xres;
803  const double yresol = yres;
804  const double zresol = zres;
805  const double pt2resol = pt2res/v.ptsq;
806 
807  mes_[label][prefix+"_ResolX"]->Fill(xresol);
808  mes_[label][prefix+"_ResolY"]->Fill(yresol);
809  mes_[label][prefix+"_ResolZ"]->Fill(zresol);
810  mes_[label][prefix+"_ResolPt2"]->Fill(pt2resol);
811 
812  mes_[label][prefix+"_ResolX_vs_PU"]->Fill(num_pileup_vertices, xresol);
813  mes_[label][prefix+"_ResolY_vs_PU"]->Fill(num_pileup_vertices, yresol);
814  mes_[label][prefix+"_ResolZ_vs_PU"]->Fill(num_pileup_vertices, zresol);
815  mes_[label][prefix+"_ResolPt2_vs_PU"]->Fill(num_pileup_vertices, pt2resol);
816 
817  mes_[label][prefix+"_ResolX_vs_NumTracks"]->Fill(v.nRecoTrk, xresol);
818  mes_[label][prefix+"_ResolY_vs_NumTracks"]->Fill(v.nRecoTrk, yresol);
819  mes_[label][prefix+"_ResolZ_vs_NumTracks"]->Fill(v.nRecoTrk, zresol);
820  mes_[label][prefix+"_ResolPt2_vs_NumTracks"]->Fill(v.nRecoTrk, pt2resol);
821 
822  mes_[label][prefix+"_ResolX_vs_Z"]->Fill(v.z, xresol);
823  mes_[label][prefix+"_ResolY_vs_Z"]->Fill(v.z, yresol);
824  mes_[label][prefix+"_ResolZ_vs_Z"]->Fill(v.z, zresol);
825  mes_[label][prefix+"_ResolPt2_vs_Z"]->Fill(v.z, pt2resol);
826 
827  mes_[label][prefix+"_ResolX_vs_Pt"]->Fill(v.pt, xresol);
828  mes_[label][prefix+"_ResolY_vs_Pt"]->Fill(v.pt, yresol);
829  mes_[label][prefix+"_ResolZ_vs_Pt"]->Fill(v.pt, zresol);
830  mes_[label][prefix+"_ResolPt2_vs_Pt"]->Fill(v.pt, pt2resol);
831 
832  mes_[label][prefix+"_PullX"]->Fill(xres/v.recVtx->xError());
833  mes_[label][prefix+"_PullY"]->Fill(yres/v.recVtx->yError());
834  mes_[label][prefix+"_PullZ"]->Fill(zres/v.recVtx->zError());
835 }
836 
838  auto found = r2s_->find(recoTrack);
839 
840  // reco track not matched to any TP
841  if(found == r2s_->end())
842  return false;
843 
844  // reco track matched to some TP from signal vertex
845  for(const auto& tp: found->val) {
846  if(tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
847  return true;
848  }
849 
850  // reco track not matched to any TP from signal vertex
851  return false;
852 }
853 
855  const std::string& label,
856  std::vector<recoPrimaryVertex>& recopvs,
857  int genpv_position_in_reco_collection,
858  bool signal_is_highest_pt) {
859  if(recopvs.empty()) return;
860 
861  std::vector<double> vtx_sumpt_sigmatched;
862  std::vector<double> vtx_sumpt2_sigmatched;
863 
864  vtx_sumpt_sigmatched.reserve(recopvs.size());
865  vtx_sumpt2_sigmatched.reserve(recopvs.size());
866 
867  // Calculate purity
868  for(auto& v: recopvs) {
869  double sumpt_all = 0;
870  double sumpt_sigmatched = 0;
871  double sumpt2_sigmatched = 0;
872  const reco::Vertex *vertex = v.recVtx;
873  for(auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
874  double pt = (*iTrack)->pt();
875  sumpt_all += pt;
876  if(matchRecoTrack2SimSignal(*iTrack)) {
877  sumpt_sigmatched += pt;
878  sumpt2_sigmatched += pt*pt;
879  }
880  }
881  v.purity = sumpt_sigmatched / sumpt_all;
882 
883  vtx_sumpt_sigmatched.push_back(sumpt_sigmatched);
884  vtx_sumpt2_sigmatched.push_back(sumpt2_sigmatched);
885  }
886 
887  const double vtxAll_sumpt_sigmatched = std::accumulate(vtx_sumpt_sigmatched.begin(), vtx_sumpt_sigmatched.end(), 0.0);
888  const double vtxNot0_sumpt_sigmatched = vtxAll_sumpt_sigmatched - vtx_sumpt_sigmatched[0];
889  const double missing = vtxNot0_sumpt_sigmatched / vtxAll_sumpt_sigmatched;
890 
891  // Fill purity
892  std::string prefix = "RecoPVAssoc2GenPVNotMatched_";
893  if(genpv_position_in_reco_collection == 0)
894  prefix = "RecoPVAssoc2GenPVMatched_";
895 
896  mes_[label][prefix+"Purity"]->Fill(recopvs[0].purity);
897  mes_[label][prefix+"Missing"]->Fill(missing);
898  auto hpurity = mes_[label][prefix+"Purity_vs_Index"];
899  for(size_t i=0; i<recopvs.size(); ++i) {
900  hpurity->Fill(i, recopvs[i].purity);
901  }
902 
903  // Fill sumpt2
904  for(size_t i=0; i<recopvs.size(); ++i) {
905  if(static_cast<int>(i) == genpv_position_in_reco_collection) {
906  mes_[label]["RecoAssoc2GenPVMatched_Pt2"]->Fill(recopvs[i].ptsq);
907  }
908  else {
909  double pt2 = recopvs[i].ptsq;
910  mes_[label]["RecoAssoc2GenPVNotMatched_Pt2"]->Fill(pt2);
911  // Subtract hard-scatter track pt2 from the pileup pt2
912  double pt2_pu = pt2-vtx_sumpt2_sigmatched[i];
913  mes_[label]["RecoAssoc2GenPVNotMatched_GenPVTracksRemoved_Pt2"]->Fill(pt2_pu);
914  }
915  }
916  if(!signal_is_highest_pt && genpv_position_in_reco_collection >= 0)
917  mes_[label]["RecoAssoc2GenPVMatchedNotHighest_Pt2"]->Fill(recopvs[genpv_position_in_reco_collection].ptsq);
918 }
919 
920 /* Extract information form TrackingParticles/TrackingVertex and fill
921  * the helper class simPrimaryVertex with proper generation-level
922  * information */
923 std::vector<PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex>
926  std::vector<PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex> simpv;
927  int current_event = -1;
928 
929  if (verbose_) {
930  std::cout << "getSimPVs TrackingVertexCollection " << std::endl;
931  }
932 
933  for (TrackingVertexCollection::const_iterator v = tVC->begin();
934  v != tVC->end(); ++v) {
935  if (verbose_) {
936  std::cout << "BunchX.EventId: " << v->eventId().bunchCrossing() << "."
937  << (v->eventId()).event() << " Position: " << v->position()
938  << " G4/HepMC Vertices: " << v->g4Vertices().size() << "/"
939  << v->genVertices().size()
940  << " t = " << v->position().t() * 1.e12
941  << " == 0:" << (v->position().t() > 0) << std::endl;
942  for (TrackingVertex::g4v_iterator gv = v->g4Vertices_begin();
943  gv != v->g4Vertices_end(); gv++) {
944  std::cout << *gv << std::endl;
945  }
946  std::cout << "----------" << std::endl;
947 
948  } // end of verbose_ session
949 
950  // I'd rather change this and select only vertices that come from
951  // BX=0. We should keep only the first vertex from all the events
952  // at BX=0.
953  if (v->eventId().bunchCrossing() != 0) continue;
954  if (v->eventId().event() != current_event) {
955  current_event = v->eventId().event();
956  } else {
957  continue;
958  }
959  // TODO(rovere) is this really necessary?
960  if (fabs(v->position().z()) > 1000) continue; // skip funny junk vertices
961 
962  // could be a new vertex, check all primaries found so far to avoid
963  // multiple entries
964  simPrimaryVertex sv(v->position().x(), v->position().y(),
965  v->position().z());
966  sv.eventId = v->eventId();
967  sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
968 
969  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin();
970  iTrack != v->daughterTracks_end(); ++iTrack) {
971  // TODO(rovere) isn't it always the case? Is it really worth
972  // checking this out?
973  // sv.eventId = (**iTrack).eventId();
974  assert((**iTrack).eventId().bunchCrossing() == 0);
975  }
976  // TODO(rovere) maybe get rid of this old logic completely ... ?
977  simPrimaryVertex* vp = NULL; // will become non-NULL if a vertex
978  // is found and then point to it
979  for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin();
980  v0 != simpv.end(); v0++) {
981  if ((sv.eventId == v0->eventId) && (fabs(sv.x - v0->x) < 1e-5) &&
982  (fabs(sv.y - v0->y) < 1e-5) && (fabs(sv.z - v0->z) < 1e-5)) {
983  vp = &(*v0);
984  break;
985  }
986  }
987  if (!vp) {
988  // this is a new vertex, add it to the list of sim-vertices
989  simpv.push_back(sv);
990  vp = &simpv.back();
991  if (verbose_) {
992  std::cout << "this is a new vertex " << sv.eventId.event() << " "
993  << sv.x << " " << sv.y << " " << sv.z << std::endl;
994  }
995  } else {
996  if (verbose_) {
997  std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " "
998  << sv.z << std::endl;
999  }
1000  }
1001 
1002  // Loop over daughter track(s) as Tracking Particles
1003  for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin();
1004  iTP != v->daughterTracks_end(); ++iTP) {
1005  auto momentum = (*(*iTP)).momentum();
1006  const reco::Track* matched_best_reco_track = nullptr;
1007  double match_quality = -1;
1008  if (use_only_charged_tracks_ && (**iTP).charge() == 0)
1009  continue;
1010  if (s2r_->find(*iTP) != s2r_->end()) {
1011  matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
1012  match_quality = (*s2r_)[*iTP][0].second;
1013  }
1014  if (verbose_) {
1015  std::cout << " Daughter momentum: " << momentum;
1016  std::cout << " Daughter type " << (*(*iTP)).pdgId();
1017  std::cout << " matched: " << (matched_best_reco_track != nullptr);
1018  std::cout << " match-quality: " << match_quality;
1019  std::cout << std::endl;
1020  }
1021  vp->ptot.setPx(vp->ptot.x() + momentum.x());
1022  vp->ptot.setPy(vp->ptot.y() + momentum.y());
1023  vp->ptot.setPz(vp->ptot.z() + momentum.z());
1024  vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
1025  vp->ptsq += ((**iTP).pt() * (**iTP).pt());
1026  // TODO(rovere) only select charged sim-particles? If so, maybe
1027  // put it as a configuration parameter?
1028  if (matched_best_reco_track) {
1030  vp->average_match_quality += match_quality;
1031  }
1032  // TODO(rovere) get rid of cuts on sim-tracks
1033  // TODO(rovere) be consistent between simulated tracks and
1034  // reconstructed tracks selection
1035  // count relevant particles
1036  if (((**iTP).pt() > 0.2) && (fabs((**iTP).eta()) < 2.5) &&
1037  (**iTP).charge() != 0) {
1038  vp->nGenTrk++;
1039  }
1040  } // End of for loop on daughters sim-particles
1041  if (vp->num_matched_reco_tracks)
1042  vp->average_match_quality /=
1043  static_cast<float>(vp->num_matched_reco_tracks);
1044  if (verbose_) {
1045  std::cout << "average number of associated tracks: "
1046  << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
1047  << " with average quality: " << vp->average_match_quality
1048  << std::endl;
1049  }
1050  } // End of for loop on tracking vertices
1051 
1052  if (verbose_) {
1053  std::cout << "------- PrimaryVertexAnalyzer4PUSlimmed simPVs from "
1054  "TrackingVertices "
1055  "-------" << std::endl;
1056  for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin();
1057  v0 != simpv.end(); v0++) {
1058  std::cout << "z=" << v0->z << " event=" << v0->eventId.event()
1059  << std::endl;
1060  }
1061  std::cout << "-----------------------------------------------" << std::endl;
1062  } // End of for summary on discovered simulated primary vertices.
1063 
1064  // In case of no simulated vertices, break here
1065  if(simpv.empty())
1066  return simpv;
1067 
1068  // Now compute the closest distance in z between all simulated vertex
1069  // first initialize
1070  auto prev_z = simpv.back().z;
1071  for(simPrimaryVertex& vsim: simpv) {
1072  vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
1073  prev_z = vsim.z;
1074  }
1075  // then calculate
1076  for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin();
1077  vsim != simpv.end(); vsim++) {
1078  std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
1079  vsim2++;
1080  for (; vsim2 != simpv.end(); vsim2++) {
1081  double distance = std::abs(vsim->z - vsim2->z);
1082  // need both to be complete
1083  vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
1084  vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
1085  }
1086  }
1087  return simpv;
1088 }
1089 
1090 /* Extract information form recoVertex and fill the helper class
1091  * recoPrimaryVertex with proper reco-level information */
1092 std::vector<PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex>
1094  const edm::Handle<edm::View<reco::Vertex>>& tVC) {
1095  std::vector<PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex> recopv;
1096 
1097  if (verbose_) {
1098  std::cout << "getRecoPVs TrackingVertexCollection " << std::endl;
1099  }
1100 
1101  for (auto v = tVC->begin(); v != tVC->end(); ++v) {
1102  if (verbose_) {
1103  std::cout << " Position: " << v->position() << std::endl;
1104  }
1105 
1106  // Skip junk vertices
1107  if (fabs(v->z()) > 1000) continue;
1108  if (v->isFake() || !v->isValid()) continue;
1109 
1110  recoPrimaryVertex sv(v->position().x(), v->position().y(),
1111  v->position().z());
1112  sv.recVtx = &(*v);
1113  sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
1114  // this is a new vertex, add it to the list of reco-vertices
1115  recopv.push_back(sv);
1117 
1118  // Loop over daughter track(s)
1119  for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
1120  auto momentum = (*(*iTrack)).innerMomentum();
1121  // TODO(rovere) better handle the pixelVertices, whose tracks
1122  // do not have the innerMomentum defined. This is a temporary
1123  // hack to overcome this problem.
1124  if (momentum.mag2() == 0)
1125  momentum = (*(*iTrack)).momentum();
1126  if (verbose_) {
1127  std::cout << " Daughter momentum: " << momentum;
1128  std::cout << std::endl;
1129  }
1130  vp->pt += std::sqrt(momentum.perp2());
1131  vp->ptsq += (momentum.perp2());
1132  vp->nRecoTrk++;
1133 
1134  auto matched = r2s_->find(*iTrack);
1135  if(matched != r2s_->end()) {
1136  vp->num_matched_sim_tracks++;
1137  }
1138 
1139  } // End of for loop on daughters reconstructed tracks
1140  } // End of for loop on tracking vertices
1141 
1142  if (verbose_) {
1143  std::cout << "------- PrimaryVertexAnalyzer4PUSlimmed recoPVs from "
1144  "VertexCollection "
1145  "-------" << std::endl;
1146  for (std::vector<recoPrimaryVertex>::iterator v0 = recopv.begin();
1147  v0 != recopv.end(); v0++) {
1148  std::cout << "z=" << v0->z << std::endl;
1149  }
1150  std::cout << "-----------------------------------------------" << std::endl;
1151  } // End of for summary on reconstructed primary vertices.
1152 
1153  // In case of no reco vertices, break here
1154  if(recopv.empty())
1155  return recopv;
1156 
1157  // Now compute the closest distance in z between all reconstructed vertex
1158  // first initialize
1159  auto prev_z = recopv.back().z;
1160  for(recoPrimaryVertex& vreco: recopv) {
1161  vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
1162  prev_z = vreco.z;
1163  }
1164  for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin();
1165  vreco != recopv.end(); vreco++) {
1166  std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
1167  vreco2++;
1168  for (; vreco2 != recopv.end(); vreco2++) {
1169  double distance = std::abs(vreco->z - vreco2->z);
1170  // need both to be complete
1171  vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
1172  vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
1173  }
1174  }
1175  return recopv;
1176 }
1177 
1179  std::vector<simPrimaryVertex> & simpv) {
1180  for (auto & v : simpv) {
1181  v.rec_vertices.clear();
1182  }
1183 }
1184 
1185 // ------------ method called to produce the data ------------
1187  std::vector<simPrimaryVertex>& simpv,
1188  const reco::VertexSimToRecoCollection& vertex_s2r) {
1189  if (verbose_) {
1190  std::cout << "PrimaryVertexAnalyzer4PUSlimmed::matchSim2RecoVertices " << std::endl;
1191  }
1192  for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin();
1193  vsim != simpv.end(); vsim++) {
1194 
1195  auto matched = vertex_s2r.find(vsim->sim_vertex);
1196  if(matched != vertex_s2r.end()) {
1197  for(const auto vertexRefQuality: matched->val) {
1198  vsim->rec_vertices.push_back(&(*(vertexRefQuality.first)));
1199  }
1200  }
1201 
1202  if (verbose_) {
1203  if (vsim->rec_vertices.size()) {
1204  for (auto const& v : vsim->rec_vertices) {
1205  std::cout << "Found a matching vertex for genVtx "
1206  << vsim->z << " at " << v->z()
1207  << " with sign: " << fabs(v->z() - vsim->z) / v->zError()
1208  << std::endl;
1209  }
1210  } else {
1211  std::cout << "No matching vertex for " << vsim->z << std::endl;
1212  }
1213  }
1214  } // end for loop on simulated vertices
1215  if (verbose_) {
1216  std::cout << "Done with matching sim vertices" << std::endl;
1217  }
1218 }
1219 
1221  std::vector<recoPrimaryVertex>& recopv,
1222  const reco::VertexRecoToSimCollection& vertex_r2s,
1223  const std::vector<simPrimaryVertex>& simpv) {
1224  for (std::vector<recoPrimaryVertex>::iterator vrec = recopv.begin();
1225  vrec != recopv.end(); vrec++) {
1226 
1227  auto matched = vertex_r2s.find(vrec->recVtxRef);
1228  if(matched != vertex_r2s.end()) {
1229  for(const auto vertexRefQuality: matched->val) {
1230  const auto tvPtr = &(*(vertexRefQuality.first));
1231  vrec->sim_vertices.push_back(tvPtr);
1232  }
1233 
1234  for(const TrackingVertex *tv: vrec->sim_vertices) {
1235  // Set pointers to internal simVertex objects
1236  for(const auto& vv: simpv) {
1237  if (&(*(vv.sim_vertex)) == tv) {
1238  vrec->sim_vertices_internal.push_back(&vv);
1239  continue;
1240  }
1241  }
1242 
1243  // Calculate number of shared tracks
1244  vrec->sim_vertices_num_shared_tracks.push_back(calculateVertexSharedTracks(*(vrec->recVtx), *tv, *r2s_));
1245  }
1246  }
1247 
1248  if (verbose_) {
1249  for (auto v : vrec->sim_vertices) {
1250  std::cout << "Found a matching vertex for reco: " << vrec->z
1251  << " at gen:" << v->position().z() << " with sign: "
1252  << fabs(vrec->z - v->position().z()) / vrec->recVtx->zError()
1253  << std::endl;
1254  }
1255  }
1256  } // end for loop on reconstructed vertices
1257 }
1258 
1260  const edm::EventSetup& iSetup) {
1261  using std::vector;
1262  using std::cout;
1263  using std::endl;
1264  using edm::Handle;
1265  using edm::View;
1266  using edm::LogInfo;
1267  using namespace reco;
1268 
1269  std::vector<float> pileUpInfo_z;
1270 
1271  // get the pileup information
1273  if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1274  for (auto const& pu_info : *puinfoH.product()) {
1275  if(do_generic_sim_plots_) {
1276  mes_["root_folder"]["GenVtx_vs_BX"]
1277  ->Fill(pu_info.getBunchCrossing(), pu_info.getPU_NumInteractions());
1278  }
1279  if (pu_info.getBunchCrossing() == 0) {
1280  pileUpInfo_z = pu_info.getPU_zpositions();
1281  if (verbose_) {
1282  for (auto const& p : pileUpInfo_z) {
1283  std::cout << "PileUpInfo on Z vertex: " << p << std::endl;
1284  }
1285  }
1286  break;
1287  }
1288  }
1289  }
1290 
1292  iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1293 
1295  iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
1296 
1297  // TODO(rovere) the idea is to put in case a track-selector in front
1298  // of this module and then use its label to get the selected tracks
1299  // out of the event instead of making an hard-coded selection in the
1300  // code.
1301 
1303  iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
1304 
1306  iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
1307 
1308  s2r_ = simToRecoH.product();
1309  r2s_ = recoToSimH.product();
1310 
1311  // Vertex associator
1313  iEvent.getByToken(vertexAssociatorToken_, vertexAssociatorH);
1314  const reco::VertexToTrackingVertexAssociator& vertexAssociator = *(vertexAssociatorH.product());
1315 
1316  std::vector<simPrimaryVertex> simpv; // a list of simulated primary
1317  // MC vertices
1318  // TODO(rovere) use move semantic?
1319  simpv = getSimPVs(TVCollectionH);
1320  // TODO(rovere) 1 vertex is not, by definition, pileup, and should
1321  // probably be subtracted?
1322  int kind_of_signal_vertex = 0;
1323  int num_pileup_vertices = simpv.size();
1324  if(do_generic_sim_plots_) mes_["root_folder"]["GenAllV_NumVertices"]->Fill(simpv.size());
1325  bool signal_is_highest_pt = std::max_element(simpv.begin(), simpv.end(),
1326  [](const simPrimaryVertex& lhs,
1327  const simPrimaryVertex& rhs) {
1328  return lhs.ptsq < rhs.ptsq;
1329  }) == simpv.begin();
1330  kind_of_signal_vertex |= (signal_is_highest_pt << HIGHEST_PT);
1331  if(do_generic_sim_plots_) {
1332  mes_["root_folder"]["SignalIsHighestPt2"]->Fill(signal_is_highest_pt ? 1. : 0.);
1333  computePairDistance(simpv, mes_["root_folder"]["GenAllV_PairDistanceZ"]);
1334  }
1335 
1336  int label_index = -1;
1337  for (auto const& vertex_token : reco_vertex_collection_tokens_) {
1338  std::vector<recoPrimaryVertex> recopv; // a list of reconstructed
1339  // primary MC vertices
1340  std::string label = reco_vertex_collections_[++label_index].label();
1342  if (!iEvent.getByToken(vertex_token, recVtxs)) {
1343  LogInfo("PrimaryVertexAnalyzer4PUSlimmed")
1344  << "Skipping vertex collection: " << label << " since it is missing."
1345  << std::endl;
1346  continue;
1347  }
1348 
1349  reco::VertexRecoToSimCollection vertex_r2s = vertexAssociator.associateRecoToSim(recVtxs, TVCollectionH);
1350  reco::VertexSimToRecoCollection vertex_s2r = vertexAssociator.associateSimToReco(recVtxs, TVCollectionH);
1351 
1352  resetSimPVAssociation(simpv);
1353  matchSim2RecoVertices(simpv, vertex_s2r);
1354  recopv = getRecoPVs(recVtxs);
1355  computePairDistance(recopv,
1356  mes_[label]["RecoAllAssoc2Gen_PairDistanceZ"]);
1357  matchReco2SimVertices(recopv, vertex_r2s, simpv);
1358 
1359  int num_total_gen_vertices_assoc2reco = 0;
1360  int num_total_reco_vertices_assoc2gen = 0;
1361  int num_total_gen_vertices_multiassoc2reco = 0;
1362  int num_total_reco_vertices_multiassoc2gen = 0;
1363  int num_total_reco_vertices_duplicate = 0;
1364  int genpv_position_in_reco_collection = -1;
1365  for (auto const& v : simpv) {
1366  float mistag = 1.;
1367  // TODO(rovere) put selectors here in front of fill* methods.
1368  if (v.eventId.event() == 0) {
1369  if (!recVtxs->empty() &&
1370  std::find(v.rec_vertices.begin(), v.rec_vertices.end(),
1371  &((*recVtxs.product())[0])) != v.rec_vertices.end()) {
1372  mistag = 0.;
1373  kind_of_signal_vertex |= (1 << IS_ASSOC2FIRST_RECO);
1374  } else {
1375  if (v.rec_vertices.size()) {
1376  kind_of_signal_vertex |= (1 << IS_ASSOC2ANY_RECO);
1377  }
1378  }
1379  mes_[label]["KindOfSignalPV"]->Fill(kind_of_signal_vertex);
1380  mes_[label]["MisTagRate"]->Fill(mistag);
1381  mes_[label]["MisTagRate_vs_PU"]->Fill(simpv.size(), mistag);
1382  mes_[label]["MisTagRate_vs_sum-pt2"]->Fill(v.ptsq, mistag);
1383  mes_[label]["MisTagRate_vs_Z"]->Fill(v.z, mistag);
1384  mes_[label]["MisTagRate_vs_R"]->Fill(v.r, mistag);
1385  mes_[label]["MisTagRate_vs_NumTracks"]->Fill(v.nGenTrk, mistag);
1386  if (signal_is_highest_pt) {
1387  mes_[label]["MisTagRateSignalIsHighest"]->Fill(mistag);
1388  mes_[label]["MisTagRateSignalIsHighest_vs_PU"]->Fill(simpv.size(),
1389  mistag);
1390  mes_[label]["MisTagRateSignalIsHighest_vs_sum-pt2"]->Fill(v.ptsq,
1391  mistag);
1392  mes_[label]["MisTagRateSignalIsHighest_vs_Z"]->Fill(v.z, mistag);
1393  mes_[label]["MisTagRateSignalIsHighest_vs_R"]->Fill(v.r, mistag);
1394  mes_[label]["MisTagRateSignalIsHighest_vs_NumTracks"]->Fill(v.nGenTrk,
1395  mistag);
1396  } else {
1397  mes_[label]["MisTagRateSignalIsNotHighest"]->Fill(mistag);
1398  mes_[label]["MisTagRateSignalIsNotHighest_vs_PU"]->Fill(simpv.size(),
1399  mistag);
1400  mes_[label]["MisTagRateSignalIsNotHighest_vs_sum-pt2"]->Fill(v.ptsq,
1401  mistag);
1402  mes_[label]["MisTagRateSignalIsNotHighest_vs_Z"]->Fill(v.z, mistag);
1403  mes_[label]["MisTagRateSignalIsNotHighest_vs_R"]->Fill(v.r, mistag);
1404  mes_[label]["MisTagRateSignalIsNotHighest_vs_NumTracks"]->
1405  Fill(v.nGenTrk, mistag);
1406  }
1407  // Now check at which location the Simulated PV has been
1408  // reconstructed in the primary vertex collection
1409  // at-hand. Mark it with fake index -1 if it was not
1410  // reconstructed at all.
1411 
1412  auto iv = (*recVtxs.product()).begin();
1413  for (int pv_position_in_reco_collection = 0;
1414  iv != (*recVtxs.product()).end();
1415  ++pv_position_in_reco_collection, ++iv) {
1416  if (std::find(v.rec_vertices.begin(), v.rec_vertices.end(),
1417  &(*iv)) != v.rec_vertices.end()) {
1418  mes_[label]["TruePVLocationIndex"]
1419  ->Fill(pv_position_in_reco_collection);
1420  const bool genPVMatchedToRecoPV = (pv_position_in_reco_collection == 0);
1421  mes_[label]["TruePVLocationIndexCumulative"]
1422  ->Fill(genPVMatchedToRecoPV ? 0 : 1);
1423 
1424  if (signal_is_highest_pt) {
1425  mes_[label]["TruePVLocationIndexSignalIsHighest"]
1426  ->Fill(pv_position_in_reco_collection);
1427  } else {
1428  mes_[label]["TruePVLocationIndexSignalIsNotHighest"]
1429  ->Fill(pv_position_in_reco_collection);
1430  }
1431 
1432  fillRecoAssociatedGenPVHistograms(label, v, genPVMatchedToRecoPV);
1433  genpv_position_in_reco_collection = pv_position_in_reco_collection;
1434  break;
1435  }
1436  }
1437 
1438  // If we reached the end, it means that the Simulated PV has not
1439  // been associated to any reconstructed vertex: mark it as
1440  // missing in the reconstructed vertex collection using the fake
1441  // index -1.
1442  if (iv == (*recVtxs.product()).end()) {
1443  mes_[label]["TruePVLocationIndex"]->Fill(-1.);
1444  mes_[label]["TruePVLocationIndexCumulative"]->Fill(-1.);
1445  if (signal_is_highest_pt)
1446  mes_[label]["TruePVLocationIndexSignalIsHighest"]->Fill(-1.);
1447  else
1448  mes_[label]["TruePVLocationIndexSignalIsNotHighest"]->Fill(-1.);
1449  }
1450  }
1451 
1452  if (v.rec_vertices.size()) num_total_gen_vertices_assoc2reco++;
1453  if (v.rec_vertices.size() > 1) num_total_gen_vertices_multiassoc2reco++;
1454  // No need to N-tplicate the Gen-related cumulative histograms:
1455  // fill them only at the first iteration
1456  if (do_generic_sim_plots_ && label_index == 0) fillGenericGenVertexHistograms(v);
1458  }
1459  calculatePurityAndFillHistograms(label, recopv, genpv_position_in_reco_collection, signal_is_highest_pt);
1460 
1461  mes_[label]["GenAllAssoc2Reco_NumVertices"]
1462  ->Fill(simpv.size(), simpv.size());
1463  mes_[label]["GenAllAssoc2RecoMatched_NumVertices"]
1464  ->Fill(simpv.size(), num_total_gen_vertices_assoc2reco);
1465  mes_[label]["GenAllAssoc2RecoMultiMatched_NumVertices"]
1466  ->Fill(simpv.size(), num_total_gen_vertices_multiassoc2reco);
1467  for (auto & v : recopv) {
1468  fillGenAssociatedRecoVertexHistograms(label, num_pileup_vertices, v);
1469  if (v.sim_vertices.size()) {
1470  num_total_reco_vertices_assoc2gen++;
1471  if (v.sim_vertices_internal[0]->rec_vertices.size() > 1) {
1472  num_total_reco_vertices_duplicate++;
1473  }
1474  }
1475  if (v.sim_vertices.size() > 1) num_total_reco_vertices_multiassoc2gen++;
1476  }
1477  mes_[label]["RecoAllAssoc2Gen_NumVertices"]
1478  ->Fill(recopv.size(), recopv.size());
1479  mes_[label]["RecoAllAssoc2GenMatched_NumVertices"]
1480  ->Fill(recopv.size(), num_total_reco_vertices_assoc2gen);
1481  mes_[label]["RecoAllAssoc2GenMultiMatched_NumVertices"]
1482  ->Fill(recopv.size(), num_total_reco_vertices_multiassoc2gen);
1483  mes_[label]["RecoAllAssoc2MultiMatchedGen_NumVertices"]
1484  ->Fill(recopv.size(), num_total_reco_vertices_duplicate);
1485  mes_[label]["RecoVtx_vs_GenVtx"]->Fill(simpv.size(), recopv.size());
1486  mes_[label]["MatchedRecoVtx_vs_GenVtx"]
1487  ->Fill(simpv.size(), num_total_reco_vertices_assoc2gen);
1488 
1489  }
1490 } // end of analyze
1491 
1492 template<class T>
1494  MonitorElement *me) {
1495  for (unsigned int i = 0; i < collection.size(); ++i) {
1496  for (unsigned int j = i+1; j < collection.size(); ++j) {
1497  me->Fill(
1498  std::abs(collection[i].z-collection[j].z));
1499  }
1500  }
1501 }
def bestMatch
Definition: deltar.py:139
T getParameter(std::string const &) const
dictionary missing
Definition: combine.py:4
void matchSim2RecoVertices(std::vector< simPrimaryVertex > &, const reco::VertexSimToRecoCollection &)
void fillGenAssociatedRecoVertexHistograms(const std::string &, int, recoPrimaryVertex &v)
std::vector< PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex > getRecoPVs(const edm::Handle< edm::View< reco::Vertex >> &)
int i
Definition: DBlmapReader.cc:9
tuple ret
prodAgent to be discontinued
int event() const
get the contents of the subdetector field (should be protected?)
bool matchRecoTrack2SimSignal(const reco::TrackBaseRef &)
void resetSimPVAssociation(std::vector< simPrimaryVertex > &)
edm::EDGetTokenT< reco::SimToRecoCollection > simToRecoAssociationToken_
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > vecPileupSummaryInfoToken_
const double xbins[]
std::vector< TrackingParticle > TrackingParticleCollection
void fillRecoAssociatedGenPVHistograms(const std::string &label, const PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex &v, bool genPVMatchedToRecoPV)
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:81
void computePairDistance(const T &collection, MonitorElement *me)
const_iterator end() const
last iterator over the map (read only)
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
double zError() const
error on z
Definition: Vertex.h:123
std::vector< SimVertex >::const_iterator g4v_iterator
assert(m_qm.get())
edm::RefToBase< reco::Vertex > VertexBaseRef
persistent reference to a Vertex, using views
Definition: VertexFwd.h:21
const_iterator find(const key_type &k) const
find element with specified reference key
virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
PrimaryVertexAnalyzer4PUSlimmed(const edm::ParameterSet &)
#define NULL
Definition: scimark2.h:8
edm::EDGetTokenT< TrackingVertexCollection > trackingVertexCollectionToken_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
void Fill(long long x)
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
int iEvent
Definition: GenABIO.cc:230
std::vector< PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex > getSimPVs(const edm::Handle< TrackingVertexCollection > &)
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
T sqrt(T t)
Definition: SSEVec.h:18
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
#define end
Definition: vmac.h:37
T min(T a, T b)
Definition: MathUtil.h:58
void calculatePurityAndFillHistograms(const std::string &, std::vector< recoPrimaryVertex > &, int, bool)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
const reco::SimToRecoCollection * s2r_
double ndof() const
Definition: Vertex.h:105
std::vector< edm::EDGetTokenT< edm::View< reco::Vertex > > > reco_vertex_collection_tokens_
void fillGenericGenVertexHistograms(const simPrimaryVertex &v)
edm::Ref< TrackingVertexCollection > TrackingVertexRef
double xError() const
error on x
Definition: Vertex.h:119
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
#define N
Definition: blowfish.cc:9
T const * product() const
Definition: Handle.h:81
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
std::vector< TrackingVertex > TrackingVertexCollection
edm::AssociationMap< edm::OneToManyWithQualityGeneric< edm::View< reco::Track >, TrackingParticleCollection, double > > RecoToSimCollection
void fillRecoAssociatedGenVertexHistograms(const std::string &, const simPrimaryVertex &v)
void matchReco2SimVertices(std::vector< recoPrimaryVertex > &, const reco::VertexRecoToSimCollection &, const std::vector< simPrimaryVertex > &)
edm::EDGetTokenT< reco::VertexToTrackingVertexAssociator > vertexAssociatorToken_
edm::EDGetTokenT< reco::RecoToSimCollection > recoToSimAssociationToken_
std::vector< edm::InputTag > reco_vertex_collections_
#define begin
Definition: vmac.h:30
std::map< std::string, std::map< std::string, MonitorElement * > > mes_
edm::AssociationMap< edm::OneToManyWithQualityGeneric< TrackingParticleCollection, edm::View< reco::Track >, double > > SimToRecoCollection
unsigned int calculateVertexSharedTracks(const reco::Vertex &recoV, const TrackingVertex &simV, const reco::RecoToSimCollection &trackRecoToSimAssociation)
tuple cout
Definition: gather_cfg.py:145
volatile std::atomic< bool > shutdown_flag false
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:76
long double T
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double yError() const
error on y
Definition: Vertex.h:121
const reco::RecoToSimCollection * r2s_
Definition: Run.h:42
edm::EDGetTokenT< TrackingParticleCollection > trackingParticleCollectionToken_
void fillResolutionAndPullHistograms(const std::string &, int, recoPrimaryVertex &v)