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