CMS 3D CMS Logo

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