CMS 3D CMS Logo

TICLCandidateValidator.cc
Go to the documentation of this file.
1 #include <numeric>
2 #include <iomanip>
3 #include <sstream>
4 
7 
9  edm::EDGetTokenT<std::vector<TICLCandidate>> ticlCandidates,
10  edm::EDGetTokenT<std::vector<TICLCandidate>> simTICLCandidatesToken,
11  edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken,
12  edm::EDGetTokenT<std::vector<ticl::Trackster>> trackstersToken,
16  bool isTICLv5)
17  : TICLCandidatesToken_(ticlCandidates),
18  simTICLCandidatesToken_(simTICLCandidatesToken),
19  recoTracksToken_(recoTracksToken),
20  trackstersToken_(trackstersToken),
21  associatorMapRtSToken_(associatorMapRtSToken),
22  associatorMapStRToken_(associatorMapStRToken),
23  associatorMapRtSPUToken_(associatorMapRtSPUToken),
24  isTICLv5_(isTICLv5) {}
25 
27 
30  std::string baseDir) const {
31  // book CAND histos
32  histograms.h_tracksters_in_candidate =
33  ibook.book1D("N of tracksters in candidate", "N of tracksters in candidate", 100, 0, 99);
34  histograms.h_candidate_raw_energy =
35  ibook.book1D("Candidates raw energy", "Candidates raw energy;E (GeV)", 250, 0, 250);
36  histograms.h_candidate_regressed_energy =
37  ibook.book1D("Candidates regressed energy", "Candidates regressed energy;E (GeV)", 250, 0, 250);
38  histograms.h_candidate_pT = ibook.book1D("Candidates pT", "Candidates pT;p_{T}", 250, 0, 250);
39  histograms.h_candidate_charge = ibook.book1D("Candidates charge", "Candidates charge;Charge", 3, -1.5, 1.5);
40  histograms.h_candidate_pdgId = ibook.book1D("Candidates PDG Id", "Candidates PDG ID", 100, -220, 220);
41  histograms.h_candidate_partType = ibook.book1D("Candidates type", "Candidates type", 9, -0.5, 8.5);
42 
43  // neutral: photon, pion, hadron
44  const std::vector<std::string> neutrals{"photons", "neutral_pions", "neutral_hadrons"};
45  for (long unsigned int i = 0; i < neutrals.size(); i++) {
46  ibook.setCurrentFolder(baseDir + "/" + neutrals[i]);
47 
48  histograms.h_neut_tracksters_in_candidate.push_back(ibook.book1D("N of tracksters in candidate for " + neutrals[i],
49  "N of tracksters in candidate for " + neutrals[i],
50  100,
51  0,
52  99));
53  histograms.h_neut_candidate_regressed_energy.push_back(ibook.book1D(
54  neutrals[i] + "candidates regressed energy", neutrals[i] + " candidates regressed energy;E (GeV)", 250, 0, 250));
55  histograms.h_neut_candidate_charge.push_back(
56  ibook.book1D(neutrals[i] + " candidates charge", neutrals[i] + " candidates charge;Charge", 3, -1.5, 1.5));
57  histograms.h_neut_candidate_pdgId.push_back(
58  ibook.book1D(neutrals[i] + " candidates PDG Id", neutrals[i] + " candidates PDG ID", 100, -220, 220));
59  histograms.h_neut_candidate_partType.push_back(
60  ibook.book1D(neutrals[i] + " candidates type", neutrals[i] + " candidates type", 9, -0.5, 8.5));
61 
62  histograms.h_den_fake_neut_energy_candidate.push_back(
63  ibook.book1D("den_fake_cand_vs_energy_" + neutrals[i], neutrals[i] + " candidates energy;E (GeV)", 50, 0, 250));
64  histograms.h_num_fake_neut_energy_candidate_pdgId.push_back(ibook.book1D(
65  "num_fake_pid_cand_vs_energy_" + neutrals[i], neutrals[i] + " PID fake vs energy;E (GeV)", 50, 0, 250));
66  histograms.h_num_fake_neut_energy_candidate_energy.push_back(
67  ibook.book1D("num_fake_energy_cand_vs_energy_" + neutrals[i],
68  neutrals[i] + " PID and energy fake vs energy;E (GeV)",
69  50,
70  0,
71  250));
72  histograms.h_den_fake_neut_pt_candidate.push_back(
73  ibook.book1D("den_fake_cand_vs_pt_" + neutrals[i], neutrals[i] + " candidates pT;p_{T} (GeV)", 50, 0, 250));
74  histograms.h_num_fake_neut_pt_candidate_pdgId.push_back(ibook.book1D(
75  "num_fake_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " PID fake vs pT;p_{T} (GeV)", 50, 0, 250));
76  histograms.h_num_fake_neut_pt_candidate_energy.push_back(
77  ibook.book1D("num_fake_energy_cand_vs_pt_" + neutrals[i],
78  neutrals[i] + " PID and energy fake vs pT;p_{T} (GeV)",
79  50,
80  0,
81  250));
82  histograms.h_den_fake_neut_eta_candidate.push_back(
83  ibook.book1D("den_fake_cand_vs_eta_" + neutrals[i], neutrals[i] + " candidates eta;#eta (GeV)", 50, -3, 3));
84  histograms.h_num_fake_neut_eta_candidate_pdgId.push_back(ibook.book1D(
85  "num_fake_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " PID fake vs eta;#eta (GeV)", 50, -3, 3));
86  histograms.h_num_fake_neut_eta_candidate_energy.push_back(
87  ibook.book1D("num_fake_energy_cand_vs_eta_" + neutrals[i],
88  neutrals[i] + " PID and energy fake vs eta;#eta (GeV)",
89  50,
90  -3,
91  3));
92  histograms.h_den_fake_neut_phi_candidate.push_back(ibook.book1D(
93  "den_fake_cand_vs_phi_" + neutrals[i], neutrals[i] + " candidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
94  histograms.h_num_fake_neut_phi_candidate_pdgId.push_back(ibook.book1D(
95  "num_fake_pid_cand_vs_phi_" + neutrals[i], neutrals[i] + " PID fake vs phi;#phi (GeV)", 50, -3.14159, 3.14159));
96  histograms.h_num_fake_neut_phi_candidate_energy.push_back(
97  ibook.book1D("num_fake_energy_cand_vs_phi_" + neutrals[i],
98  neutrals[i] + " PID and energy fake vs phi;#phi (GeV)",
99  50,
100  -3.14159,
101  3.14159));
102 
103  histograms.h_den_neut_energy_candidate.push_back(
104  ibook.book1D("den_cand_vs_energy_" + neutrals[i], neutrals[i] + " simCandidates energy;E (GeV)", 50, 0, 250));
105  histograms.h_num_neut_energy_candidate_pdgId.push_back(
106  ibook.book1D("num_pid_cand_vs_energy_" + neutrals[i],
107  neutrals[i] + " track and PID efficiency vs energy;E (GeV)",
108  50,
109  0,
110  250));
111  histograms.h_num_neut_energy_candidate_energy.push_back(
112  ibook.book1D("num_energy_cand_vs_energy_" + neutrals[i],
113  neutrals[i] + " track, PID and energy efficiency vs energy;E (GeV)",
114  50,
115  0,
116  250));
117  histograms.h_den_neut_pt_candidate.push_back(
118  ibook.book1D("den_cand_vs_pt_" + neutrals[i], neutrals[i] + " simCandidates pT;p_{T} (GeV)", 50, 0, 250));
119  histograms.h_num_neut_pt_candidate_pdgId.push_back(ibook.book1D(
120  "num_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
121  histograms.h_num_neut_pt_candidate_energy.push_back(
122  ibook.book1D("num_energy_cand_vs_pt_" + neutrals[i],
123  neutrals[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)",
124  50,
125  0,
126  250));
127  histograms.h_den_neut_eta_candidate.push_back(
128  ibook.book1D("den_cand_vs_eta_" + neutrals[i], neutrals[i] + " simCandidates eta;#eta (GeV)", 50, -3, 3));
129  histograms.h_num_neut_eta_candidate_pdgId.push_back(ibook.book1D(
130  "num_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " track and PID efficiency vs eta;#eta (GeV)", 50, -3, 3));
131  histograms.h_num_neut_eta_candidate_energy.push_back(
132  ibook.book1D("num_energy_cand_vs_eta_" + neutrals[i],
133  neutrals[i] + " track, PID and energy efficiency vs eta;#eta (GeV)",
134  50,
135  -3,
136  3));
137  histograms.h_den_neut_phi_candidate.push_back(ibook.book1D(
138  "den_cand_vs_phi_" + neutrals[i], neutrals[i] + " simCandidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
139  histograms.h_num_neut_phi_candidate_pdgId.push_back(
140  ibook.book1D("num_pid_cand_vs_phi_" + neutrals[i],
141  neutrals[i] + " track and PID efficiency vs phi;#phi (GeV)",
142  50,
143  -3.14159,
144  3.14159));
145  histograms.h_num_neut_phi_candidate_energy.push_back(
146  ibook.book1D("num_energy_cand_vs_phi_" + neutrals[i],
147  neutrals[i] + " track, PID and energy efficiency vs phi;#phi (GeV)",
148  50,
149  -3.14159,
150  3.14159));
151  }
152  // charged: electron, muon, hadron
153  const std::vector<std::string> charged{"electrons", "muons", "charged_hadrons"};
154  for (long unsigned int i = 0; i < charged.size(); i++) {
155  ibook.setCurrentFolder(baseDir + "/" + charged[i]);
156 
157  histograms.h_chg_tracksters_in_candidate.push_back(ibook.book1D(
158  "N of tracksters in candidate for " + charged[i], "N of tracksters in candidate for " + charged[i], 100, 0, 99));
159  histograms.h_chg_candidate_regressed_energy.push_back(ibook.book1D(
160  charged[i] + "candidates regressed energy", charged[i] + " candidates regressed energy;E (GeV)", 250, 0, 250));
161  histograms.h_chg_candidate_charge.push_back(
162  ibook.book1D(charged[i] + " candidates charge", charged[i] + " candidates charge;Charge", 3, -1.5, 1.5));
163  histograms.h_chg_candidate_pdgId.push_back(
164  ibook.book1D(charged[i] + " candidates PDG Id", charged[i] + " candidates PDG ID", 100, -220, 220));
165  histograms.h_chg_candidate_partType.push_back(
166  ibook.book1D(charged[i] + " candidates type", charged[i] + " candidates type", 9, -0.5, 8.5));
167 
168  histograms.h_den_fake_chg_energy_candidate.push_back(
169  ibook.book1D("den_fake_cand_vs_energy_" + charged[i], charged[i] + " candidates energy;E (GeV)", 50, 0, 250));
170  histograms.h_num_fake_chg_energy_candidate_track.push_back(ibook.book1D(
171  "num_fake_track_cand_vs_energy_" + charged[i], charged[i] + " track fake vs energy;E (GeV)", 50, 0, 250));
172  histograms.h_num_fake_chg_energy_candidate_pdgId.push_back(ibook.book1D(
173  "num_fake_pid_cand_vs_energy_" + charged[i], charged[i] + " track and PID fake vs energy;E (GeV)", 50, 0, 250));
174  histograms.h_num_fake_chg_energy_candidate_energy.push_back(
175  ibook.book1D("num_fake_energy_cand_vs_energy_" + charged[i],
176  charged[i] + " track, PID and energy fake vs energy;E (GeV)",
177  50,
178  0,
179  250));
180  histograms.h_den_fake_chg_pt_candidate.push_back(
181  ibook.book1D("den_fake_cand_vs_pt_" + charged[i], charged[i] + " candidates pT;p_{T} (GeV)", 50, 0, 250));
182  histograms.h_num_fake_chg_pt_candidate_track.push_back(ibook.book1D(
183  "num_fake_track_cand_vs_pt_" + charged[i], charged[i] + " track fake vs pT;p_{T} (GeV)", 50, 0, 250));
184  histograms.h_num_fake_chg_pt_candidate_pdgId.push_back(ibook.book1D(
185  "num_fake_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID fake vs pT;p_{T} (GeV)", 50, 0, 250));
186  histograms.h_num_fake_chg_pt_candidate_energy.push_back(
187  ibook.book1D("num_fake_energy_cand_vs_pt_" + charged[i],
188  charged[i] + " track, PID and energy fake vs pT;p_{T} (GeV)",
189  50,
190  0,
191  250));
192  histograms.h_den_fake_chg_eta_candidate.push_back(
193  ibook.book1D("den_fake_cand_vs_eta_" + charged[i], charged[i] + " candidates eta;#eta (GeV)", 50, -3, 3));
194  histograms.h_num_fake_chg_eta_candidate_track.push_back(ibook.book1D(
195  "num_fake_track_cand_vs_eta_" + charged[i], charged[i] + " track fake vs eta;#eta (GeV)", 50, -3, 3));
196  histograms.h_num_fake_chg_eta_candidate_pdgId.push_back(ibook.book1D(
197  "num_fake_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID fake vs eta;#eta (GeV)", 50, -3, 3));
198  histograms.h_num_fake_chg_eta_candidate_energy.push_back(
199  ibook.book1D("num_fake_energy_cand_vs_eta_" + charged[i],
200  charged[i] + " track, PID and energy fake vs eta;#eta (GeV)",
201  50,
202  -3,
203  3));
204  histograms.h_den_fake_chg_phi_candidate.push_back(ibook.book1D(
205  "den_fake_cand_vs_phi_" + charged[i], charged[i] + " candidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
206  histograms.h_num_fake_chg_phi_candidate_track.push_back(ibook.book1D("num_fake_track_cand_vs_phi_" + charged[i],
207  charged[i] + " track fake vs phi;#phi (GeV)",
208  50,
209  -3.14159,
210  3.14159));
211  histograms.h_num_fake_chg_phi_candidate_pdgId.push_back(
212  ibook.book1D("num_fake_pid_cand_vs_phi_" + charged[i],
213  charged[i] + " track and PID fake vs phi;#phi (GeV)",
214  50,
215  -3.14159,
216  3.14159));
217  histograms.h_num_fake_chg_phi_candidate_energy.push_back(
218  ibook.book1D("num_fake_energy_cand_vs_phi_" + charged[i],
219  charged[i] + " track, PID and energy fake vs phi;#phi (GeV)",
220  50,
221  -3.14159,
222  3.14159));
223 
224  histograms.h_den_chg_energy_candidate.push_back(
225  ibook.book1D("den_cand_vs_energy_" + charged[i], charged[i] + " simCandidates energy;E (GeV)", 50, 0, 250));
226  histograms.h_num_chg_energy_candidate_track.push_back(ibook.book1D(
227  "num_track_cand_vs_energy_" + charged[i], charged[i] + " track efficiency vs energy;E (GeV)", 50, 0, 250));
228  histograms.h_num_chg_energy_candidate_pdgId.push_back(ibook.book1D(
229  "num_pid_cand_vs_energy_" + charged[i], charged[i] + " track and PID efficiency vs energy;E (GeV)", 50, 0, 250));
230  histograms.h_num_chg_energy_candidate_energy.push_back(
231  ibook.book1D("num_energy_cand_vs_energy_" + charged[i],
232  charged[i] + " track, PID and energy efficiency vs energy;E (GeV)",
233  50,
234  0,
235  250));
236  histograms.h_den_chg_pt_candidate.push_back(
237  ibook.book1D("den_cand_vs_pt_" + charged[i], charged[i] + " simCandidates pT;p_{T} (GeV)", 50, 0, 250));
238  histograms.h_num_chg_pt_candidate_track.push_back(ibook.book1D(
239  "num_track_cand_vs_pt_" + charged[i], charged[i] + " track efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
240  histograms.h_num_chg_pt_candidate_pdgId.push_back(ibook.book1D(
241  "num_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
242  histograms.h_num_chg_pt_candidate_energy.push_back(
243  ibook.book1D("num_energy_cand_vs_pt_" + charged[i],
244  charged[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)",
245  50,
246  0,
247  250));
248  histograms.h_den_chg_eta_candidate.push_back(
249  ibook.book1D("den_cand_vs_eta_" + charged[i], charged[i] + " simCandidates eta;#eta (GeV)", 50, -3, 3));
250  histograms.h_num_chg_eta_candidate_track.push_back(ibook.book1D(
251  "num_track_cand_vs_eta_" + charged[i], charged[i] + " track efficiency vs eta;#eta (GeV)", 50, -3, 3));
252  histograms.h_num_chg_eta_candidate_pdgId.push_back(ibook.book1D(
253  "num_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID efficiency vs eta;#eta (GeV)", 50, -3, 3));
254  histograms.h_num_chg_eta_candidate_energy.push_back(
255  ibook.book1D("num_energy_cand_vs_eta_" + charged[i],
256  charged[i] + " track, PID and energy efficiency vs eta;#eta (GeV)",
257  50,
258  -3,
259  3));
260  histograms.h_den_chg_phi_candidate.push_back(ibook.book1D(
261  "den_cand_vs_phi_" + charged[i], charged[i] + " simCandidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
262  histograms.h_num_chg_phi_candidate_track.push_back(ibook.book1D("num_track_cand_vs_phi_" + charged[i],
263  charged[i] + " track efficiency vs phi;#phi (GeV)",
264  50,
265  -3.14159,
266  3.14159));
267  histograms.h_num_chg_phi_candidate_pdgId.push_back(
268  ibook.book1D("num_pid_cand_vs_phi_" + charged[i],
269  charged[i] + " track and PID efficiency vs phi;#phi (GeV)",
270  50,
271  -3.14159,
272  3.14159));
273  histograms.h_num_chg_phi_candidate_energy.push_back(
274  ibook.book1D("num_energy_cand_vs_phi_" + charged[i],
275  charged[i] + " track, PID and energy efficiency vs phi;#phi (GeV)",
276  50,
277  -3.14159,
278  3.14159));
279  }
280 }
281 
283  const Histograms& histograms,
284  edm::Handle<ticl::TracksterCollection> simTrackstersCP_h) const {
285  auto TICLCandidates = event.get(TICLCandidatesToken_);
286 
287  edm::Handle<std::vector<TICLCandidate>> simTICLCandidates_h;
288  event.getByToken(simTICLCandidatesToken_, simTICLCandidates_h);
289  auto simTICLCandidates = *simTICLCandidates_h;
290 
292  event.getByToken(recoTracksToken_, recoTracks_h);
293  auto recoTracks = *recoTracks_h;
294 
296  event.getByToken(trackstersToken_, Tracksters_h);
297  auto trackstersMerged = *Tracksters_h;
298 
300  event.getByToken(associatorMapRtSToken_, mergeTsRecoToSim_h);
301  auto const& mergeTsRecoToSimMap = *mergeTsRecoToSim_h;
302 
304  event.getByToken(associatorMapStRToken_, mergeTsSimToReco_h);
305  auto const& mergeTsSimToRecoMap = *mergeTsSimToReco_h;
306 
308  event.getByToken(associatorMapRtSPUToken_, mergeTsRecoToSimPU_h);
309  auto const& mergeTsRecoToSimPUMap = *mergeTsRecoToSimPU_h;
310 
311  // candidates plots
312  for (const auto& cand : TICLCandidates) {
313  histograms.h_tracksters_in_candidate->Fill(cand.tracksters().size());
314  histograms.h_candidate_raw_energy->Fill(cand.rawEnergy());
315  histograms.h_candidate_regressed_energy->Fill(cand.energy());
316  histograms.h_candidate_pT->Fill(cand.pt());
317  histograms.h_candidate_charge->Fill(cand.charge());
318  histograms.h_candidate_pdgId->Fill(cand.pdgId());
319  const auto& arr = cand.idProbabilities();
320  histograms.h_candidate_partType->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
321  }
322 
323  std::vector<int> chargedCandidates;
324  std::vector<int> neutralCandidates;
325  chargedCandidates.reserve(simTICLCandidates.size());
326  neutralCandidates.reserve(simTICLCandidates.size());
327 
328  for (size_t i = 0; i < simTICLCandidates.size(); ++i) {
329  const auto& simCand = simTICLCandidates[i];
330  const auto particleType = ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), simCand.charge());
334  chargedCandidates.emplace_back(i);
338  neutralCandidates.emplace_back(i);
339  // should consider also unknown ?
340  }
341 
342  chargedCandidates.shrink_to_fit();
343  neutralCandidates.shrink_to_fit();
344 
345  for (const auto i : chargedCandidates) {
346  const auto& simCand = simTICLCandidates[i];
347  auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 1)));
348  /* 11 (type 1) becomes 0
349  * 13 (type 2) becomes 1
350  * 211 (type 4) becomes 2
351  */
352  int32_t simCandTrackIdx = -1;
353  if (simCand.trackPtr().get() != nullptr)
354  simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
355  else {
356  // no reco track, but simCand is charged
357  continue;
358  }
359  if (simCand.trackPtr().get()->pt() < 1 or simCand.trackPtr().get()->missingOuterHits() > 5 or
360  not simCand.trackPtr().get()->quality(reco::TrackBase::highPurity))
361  continue;
362 
363  // +1 to all denominators
364  histograms.h_den_chg_energy_candidate[index]->Fill(simCand.rawEnergy());
365  histograms.h_den_chg_pt_candidate[index]->Fill(simCand.pt());
366  histograms.h_den_chg_eta_candidate[index]->Fill(simCand.eta());
367  histograms.h_den_chg_phi_candidate[index]->Fill(simCand.phi());
368 
369  int32_t cand_idx = -1;
370  const edm::Ref<ticl::TracksterCollection> stsRef(simTrackstersCP_h, i);
371  const auto ts_iter = mergeTsSimToRecoMap.find(stsRef);
372  float shared_energy = 0.;
373  // search for reco cand associated
374  if (ts_iter != mergeTsSimToRecoMap.end()) {
375  const auto& tsAssoc = (ts_iter->val);
376  std::vector<uint32_t> MergeTracksters_simToReco;
377  std::vector<float> MergeTracksters_simToReco_score;
378  std::vector<float> MergeTracksters_simToReco_sharedE;
379  MergeTracksters_simToReco.reserve(tsAssoc.size());
380  MergeTracksters_simToReco_score.reserve(tsAssoc.size());
381  MergeTracksters_simToReco_sharedE.reserve(tsAssoc.size());
382  for (auto& ts : tsAssoc) {
383  auto ts_id = (ts.first).get() - (edm::Ref<ticl::TracksterCollection>(Tracksters_h, 0)).get();
384  MergeTracksters_simToReco.push_back(ts_id);
385  MergeTracksters_simToReco_score.push_back(ts.second.second);
386  MergeTracksters_simToReco_sharedE.push_back(ts.second.first);
387  }
388  auto min_idx = std::min_element(MergeTracksters_simToReco_score.begin(), MergeTracksters_simToReco_score.end());
389  if (*min_idx != 1) {
390  cand_idx = MergeTracksters_simToReco[min_idx - MergeTracksters_simToReco_score.begin()];
391  shared_energy = MergeTracksters_simToReco_sharedE[min_idx - MergeTracksters_simToReco_score.begin()];
392  }
393  }
394 
395  // no reco associated to sim
396  if (cand_idx == -1)
397  continue;
398 
399  auto& recoCand = TICLCandidates[cand_idx];
400  if (isTICLv5_) {
401  // cand_idx is the tsMerge index, find the ts in the candidates collection
402  auto const tsPtr = edm::Ptr<ticl::Trackster>(Tracksters_h, cand_idx);
403  auto cand_it = std::find_if(TICLCandidates.begin(), TICLCandidates.end(), [tsPtr](TICLCandidate const& cand) {
404  if (!cand.tracksters().empty())
405  return cand.tracksters()[0] == tsPtr;
406  else
407  return false;
408  });
409  if (cand_it != TICLCandidates.end())
410  recoCand = *cand_it;
411  else
412  continue;
413  }
414 
415  if (recoCand.trackPtr().get() != nullptr) {
416  const auto candTrackIdx = recoCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
417  if (simCandTrackIdx == candTrackIdx) {
418  // +1 to track num
419  histograms.h_num_chg_energy_candidate_track[index]->Fill(simCand.rawEnergy());
420  histograms.h_num_chg_pt_candidate_track[index]->Fill(simCand.pt());
421  histograms.h_num_chg_eta_candidate_track[index]->Fill(simCand.eta());
422  histograms.h_num_chg_phi_candidate_track[index]->Fill(simCand.phi());
423  } else {
424  continue;
425  }
426  } else {
427  continue;
428  }
429 
430  //step 2: PID
431  if (simCand.pdgId() == recoCand.pdgId()) {
432  // +1 to num pdg id
433  histograms.h_num_chg_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy());
434  histograms.h_num_chg_pt_candidate_pdgId[index]->Fill(simCand.pt());
435  histograms.h_num_chg_eta_candidate_pdgId[index]->Fill(simCand.eta());
436  histograms.h_num_chg_phi_candidate_pdgId[index]->Fill(simCand.phi());
437 
438  //step 3: energy
439  if (shared_energy / simCand.rawEnergy() > 0.5) {
440  // +1 to ene num
441  histograms.h_num_chg_energy_candidate_energy[index]->Fill(simCand.rawEnergy());
442  histograms.h_num_chg_pt_candidate_energy[index]->Fill(simCand.pt());
443  histograms.h_num_chg_eta_candidate_energy[index]->Fill(simCand.eta());
444  histograms.h_num_chg_phi_candidate_energy[index]->Fill(simCand.phi());
445  }
446  }
447  }
448 
449  for (const auto i : neutralCandidates) {
450  const auto& simCand = simTICLCandidates[i];
451  auto index = int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 0)) / 2;
452  /* 22 (type 0) becomes 0
453  * 111 (type 3) becomes 1
454  * 130 (type 5) becomes 2
455  */
456  histograms.h_den_neut_energy_candidate[index]->Fill(simCand.rawEnergy());
457  histograms.h_den_neut_pt_candidate[index]->Fill(simCand.pt());
458  histograms.h_den_neut_eta_candidate[index]->Fill(simCand.eta());
459  histograms.h_den_neut_phi_candidate[index]->Fill(simCand.phi());
460 
461  int32_t cand_idx = -1;
462  const edm::Ref<ticl::TracksterCollection> stsRef(simTrackstersCP_h, i);
463  const auto ts_iter = mergeTsSimToRecoMap.find(stsRef);
464  float shared_energy = 0.;
465  // search for reco cand associated
466  if (ts_iter != mergeTsSimToRecoMap.end()) {
467  const auto& tsAssoc = (ts_iter->val);
468  std::vector<uint32_t> MergeTracksters_simToReco;
469  std::vector<float> MergeTracksters_simToReco_score;
470  std::vector<float> MergeTracksters_simToReco_sharedE;
471  MergeTracksters_simToReco.reserve(tsAssoc.size());
472  MergeTracksters_simToReco_score.reserve(tsAssoc.size());
473  MergeTracksters_simToReco_sharedE.reserve(tsAssoc.size());
474  for (auto& ts : tsAssoc) {
475  auto ts_id = (ts.first).get() - (edm::Ref<ticl::TracksterCollection>(Tracksters_h, 0)).get();
476  MergeTracksters_simToReco.push_back(ts_id);
477  MergeTracksters_simToReco_score.push_back(ts.second.second);
478  MergeTracksters_simToReco_sharedE.push_back(ts.second.first);
479  }
480  auto min_idx = std::min_element(MergeTracksters_simToReco_score.begin(), MergeTracksters_simToReco_score.end());
481  if (*min_idx != 1) {
482  cand_idx = MergeTracksters_simToReco[min_idx - MergeTracksters_simToReco_score.begin()];
483  shared_energy = MergeTracksters_simToReco_sharedE[min_idx - MergeTracksters_simToReco_score.begin()];
484  }
485  }
486 
487  // no reco associated to sim
488  if (cand_idx == -1)
489  continue;
490 
491  auto& recoCand = TICLCandidates[cand_idx];
492  if (isTICLv5_) {
493  // cand_idx is the tsMerge index, find the ts in the candidates collection
494  auto const tsPtr = edm::Ptr<ticl::Trackster>(Tracksters_h, cand_idx);
495  auto cand_it = std::find_if(TICLCandidates.begin(), TICLCandidates.end(), [tsPtr](TICLCandidate const& cand) {
496  if (!cand.tracksters().empty())
497  return cand.tracksters()[0] == tsPtr;
498  else
499  return false;
500  });
501  if (cand_it != TICLCandidates.end())
502  recoCand = *cand_it;
503  else
504  continue;
505  }
506 
507  if (recoCand.trackPtr().get() != nullptr)
508  continue;
509 
510  //step 2: PID
511  if (simCand.pdgId() == recoCand.pdgId()) {
512  // +1 to num pdg id
513  histograms.h_num_neut_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy());
514  histograms.h_num_neut_pt_candidate_pdgId[index]->Fill(simCand.pt());
515  histograms.h_num_neut_eta_candidate_pdgId[index]->Fill(simCand.eta());
516  histograms.h_num_neut_phi_candidate_pdgId[index]->Fill(simCand.phi());
517 
518  //step 3: energy
519  if (shared_energy / simCand.rawEnergy() > 0.5) {
520  // +1 to ene num
521  histograms.h_num_neut_energy_candidate_energy[index]->Fill(simCand.rawEnergy());
522  histograms.h_num_neut_pt_candidate_energy[index]->Fill(simCand.pt());
523  histograms.h_num_neut_eta_candidate_energy[index]->Fill(simCand.eta());
524  histograms.h_num_neut_phi_candidate_energy[index]->Fill(simCand.phi());
525  }
526  }
527  }
528 
529  // FAKE rate
530  chargedCandidates.clear();
531  neutralCandidates.clear();
532  chargedCandidates.reserve(TICLCandidates.size());
533  neutralCandidates.reserve(TICLCandidates.size());
534 
535  auto isCharged = [](int pdgId) {
536  pdgId = std::abs(pdgId);
537  return (pdgId == 11 or pdgId == 211 or pdgId == 13);
538  };
539 
540  for (size_t i = 0; i < TICLCandidates.size(); ++i) {
541  const auto& cand = TICLCandidates[i];
542  const auto& charged = isCharged(cand.pdgId());
543  if (charged)
544  chargedCandidates.emplace_back(i);
545  else
546  neutralCandidates.emplace_back(i);
547 
548  // should consider also unknown ?
549  }
550 
551  chargedCandidates.shrink_to_fit();
552  neutralCandidates.shrink_to_fit();
553 
554  // loop on charged
555  for (const auto i : chargedCandidates) {
556  const auto& cand = TICLCandidates[i];
557  auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 1)));
558  /* 11 (type 1) becomes 0
559  * 13 (type 2) becomes 1
560  * 211 (type 4) becomes 2
561  */
562  int32_t candTrackIdx = -1;
563  candTrackIdx = cand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
564 
565  if (cand.tracksters().empty())
566  continue;
567 
568  // i is the candidate idx == ts idx only in v4, find ts_idx in v5
569  auto mergeTs_id = i;
570  if (isTICLv5_) {
571  mergeTs_id = cand.tracksters()[0].get() - edm::Ptr<ticl::Trackster>(Tracksters_h, 0).get();
572  }
573  // remove PU tracksters
574  const edm::Ref<ticl::TracksterCollection> tsRef(Tracksters_h, mergeTs_id);
575  auto const sts_iterPU = mergeTsRecoToSimPUMap.find(tsRef);
576  if (sts_iterPU != mergeTsRecoToSimPUMap.end()) {
577  const auto& stsPUAssociated = sts_iterPU->val;
578  if (stsPUAssociated[0].second.first / (*Tracksters_h)[mergeTs_id].raw_energy() > 0.95)
579  continue;
580  }
581 
582  // +1 to all denominators
583  histograms.h_den_fake_chg_energy_candidate[index]->Fill(cand.rawEnergy());
584  histograms.h_den_fake_chg_pt_candidate[index]->Fill(cand.pt());
585  histograms.h_den_fake_chg_eta_candidate[index]->Fill(cand.eta());
586  histograms.h_den_fake_chg_phi_candidate[index]->Fill(cand.phi());
587 
588  histograms.h_chg_tracksters_in_candidate[index]->Fill(cand.tracksters().size());
589  histograms.h_chg_candidate_regressed_energy[index]->Fill(cand.energy());
590  histograms.h_chg_candidate_charge[index]->Fill(cand.charge());
591  histograms.h_chg_candidate_pdgId[index]->Fill(cand.pdgId());
592  const auto& arr = cand.idProbabilities();
593  histograms.h_chg_candidate_partType[index]->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
594 
595  int32_t simCand_idx = -1;
596  const auto sts_iter = mergeTsRecoToSimMap.find(tsRef);
597  float shared_energy = 0.;
598  // search for reco cand associated
599  if (sts_iter != mergeTsRecoToSimMap.end()) {
600  const auto& stsAssoc = (sts_iter->val);
601  std::vector<uint32_t> MergeTracksters_recoToSim;
602  std::vector<float> MergeTracksters_recoToSim_score;
603  std::vector<float> MergeTracksters_recoToSim_sharedE;
604  MergeTracksters_recoToSim.reserve(stsAssoc.size());
605  MergeTracksters_recoToSim_score.reserve(stsAssoc.size());
606  MergeTracksters_recoToSim_sharedE.reserve(stsAssoc.size());
607  for (auto& sts : stsAssoc) {
608  auto sts_id = (sts.first).get() - (edm::Ref<ticl::TracksterCollection>(simTrackstersCP_h, 0)).get();
609  MergeTracksters_recoToSim.push_back(sts_id);
610  MergeTracksters_recoToSim_score.push_back(sts.second.second);
611  MergeTracksters_recoToSim_sharedE.push_back(sts.second.first);
612  }
613  auto min_idx = std::min_element(MergeTracksters_recoToSim_score.begin(), MergeTracksters_recoToSim_score.end());
614  if (*min_idx != 1) {
615  simCand_idx = MergeTracksters_recoToSim[min_idx - MergeTracksters_recoToSim_score.begin()];
616  shared_energy = MergeTracksters_recoToSim_sharedE[min_idx - MergeTracksters_recoToSim_score.begin()];
617  }
618  }
619 
620  if (simCand_idx == -1)
621  continue;
622 
623  const auto& simCand = simTICLCandidates[simCand_idx];
624  if (simCand.trackPtr().get() != nullptr) {
625  const auto simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
626  if (simCandTrackIdx != candTrackIdx) {
627  // fake += 1
628  histograms.h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy());
629  histograms.h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt());
630  histograms.h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta());
631  histograms.h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi());
632  continue;
633  }
634  } else {
635  // fake += 1
636  histograms.h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy());
637  histograms.h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt());
638  histograms.h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta());
639  histograms.h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi());
640  continue;
641  }
642 
643  //step 2: PID
644  if (simCand.pdgId() != cand.pdgId()) {
645  // +1 to num fake pdg id
646  histograms.h_num_fake_chg_energy_candidate_pdgId[index]->Fill(cand.rawEnergy());
647  histograms.h_num_fake_chg_pt_candidate_pdgId[index]->Fill(cand.pt());
648  histograms.h_num_fake_chg_eta_candidate_pdgId[index]->Fill(cand.eta());
649  histograms.h_num_fake_chg_phi_candidate_pdgId[index]->Fill(cand.phi());
650  continue;
651  }
652 
653  //step 3: energy
654  if (shared_energy / simCand.rawEnergy() < 0.5) {
655  // +1 to ene num
656  histograms.h_num_fake_chg_energy_candidate_energy[index]->Fill(cand.rawEnergy());
657  histograms.h_num_fake_chg_pt_candidate_energy[index]->Fill(cand.pt());
658  histograms.h_num_fake_chg_eta_candidate_energy[index]->Fill(cand.eta());
659  histograms.h_num_fake_chg_phi_candidate_energy[index]->Fill(cand.phi());
660  }
661  }
662  // loop on neutrals
663  for (const auto i : neutralCandidates) {
664  const auto& cand = TICLCandidates[i];
665  auto index = int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 0)) / 2;
666  /* 22 (type 0) becomes 0
667  * 111 (type 3) becomes 1
668  * 130 (type 5) becomes 2
669  */
670 
671  if (cand.tracksters().empty())
672  continue;
673 
674  // i is the candidate idx == ts idx only in v4, find ts_idx in v5
675  auto mergeTs_id = i;
676  if (isTICLv5_) {
677  mergeTs_id = cand.tracksters()[0].get() - edm::Ptr<ticl::Trackster>(Tracksters_h, 0).get();
678  }
679  // remove PU tracksters
680  const edm::Ref<ticl::TracksterCollection> tsRef(Tracksters_h, mergeTs_id);
681  auto const sts_iterPU = mergeTsRecoToSimPUMap.find(tsRef);
682  if (sts_iterPU != mergeTsRecoToSimPUMap.end()) {
683  const auto& stsPUAssociated = sts_iterPU->val;
684  if (stsPUAssociated[0].second.first / (*Tracksters_h)[mergeTs_id].raw_energy() > 0.95)
685  continue;
686  }
687 
688  // +1 to all denominators
689  histograms.h_den_fake_neut_energy_candidate[index]->Fill(cand.rawEnergy());
690  histograms.h_den_fake_neut_pt_candidate[index]->Fill(cand.pt());
691  histograms.h_den_fake_neut_eta_candidate[index]->Fill(cand.eta());
692  histograms.h_den_fake_neut_phi_candidate[index]->Fill(cand.phi());
693 
694  histograms.h_neut_tracksters_in_candidate[index]->Fill(cand.tracksters().size());
695  histograms.h_neut_candidate_regressed_energy[index]->Fill(cand.energy());
696  histograms.h_neut_candidate_charge[index]->Fill(cand.charge());
697  histograms.h_neut_candidate_pdgId[index]->Fill(cand.pdgId());
698  const auto& arr = cand.idProbabilities();
699  histograms.h_neut_candidate_partType[index]->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
700 
701  int32_t simCand_idx = -1;
702  const auto sts_iter = mergeTsRecoToSimMap.find(tsRef);
703  float shared_energy = 0.;
704  // search for reco cand associated
705  if (sts_iter != mergeTsRecoToSimMap.end()) {
706  const auto& stsAssoc = (sts_iter->val);
707  std::vector<uint32_t> MergeTracksters_recoToSim;
708  std::vector<float> MergeTracksters_recoToSim_score;
709  std::vector<float> MergeTracksters_recoToSim_sharedE;
710  MergeTracksters_recoToSim.reserve(stsAssoc.size());
711  MergeTracksters_recoToSim_score.reserve(stsAssoc.size());
712  MergeTracksters_recoToSim_sharedE.reserve(stsAssoc.size());
713  for (auto& sts : stsAssoc) {
714  auto sts_id = (sts.first).get() - (edm::Ref<ticl::TracksterCollection>(simTrackstersCP_h, 0)).get();
715  MergeTracksters_recoToSim.push_back(sts_id);
716  MergeTracksters_recoToSim_score.push_back(sts.second.second);
717  MergeTracksters_recoToSim_sharedE.push_back(sts.second.first);
718  }
719  auto min_idx = std::min_element(MergeTracksters_recoToSim_score.begin(), MergeTracksters_recoToSim_score.end());
720  if (*min_idx != 1) {
721  simCand_idx = MergeTracksters_recoToSim[min_idx - MergeTracksters_recoToSim_score.begin()];
722  shared_energy = MergeTracksters_recoToSim_sharedE[min_idx - MergeTracksters_recoToSim_score.begin()];
723  }
724  }
725 
726  if (simCand_idx == -1)
727  continue;
728 
729  const auto& simCand = simTICLCandidates[simCand_idx];
730 
731  //step 2: PID
732  if (simCand.pdgId() != cand.pdgId()) {
733  // +1 to num fake pdg id
734  histograms.h_num_fake_neut_energy_candidate_pdgId[index]->Fill(cand.rawEnergy());
735  histograms.h_num_fake_neut_pt_candidate_pdgId[index]->Fill(cand.pt());
736  histograms.h_num_fake_neut_eta_candidate_pdgId[index]->Fill(cand.eta());
737  histograms.h_num_fake_neut_phi_candidate_pdgId[index]->Fill(cand.phi());
738  continue;
739  }
740 
741  //step 3: energy
742  if (shared_energy / simCand.rawEnergy() < 0.5) {
743  // +1 to ene num
744  histograms.h_num_fake_neut_energy_candidate_energy[index]->Fill(cand.rawEnergy());
745  histograms.h_num_fake_neut_pt_candidate_energy[index]->Fill(cand.pt());
746  histograms.h_num_fake_neut_eta_candidate_energy[index]->Fill(cand.eta());
747  histograms.h_num_fake_neut_phi_candidate_energy[index]->Fill(cand.phi());
748  }
749  }
750 }
edm::EDGetTokenT< std::vector< TICLCandidate > > simTICLCandidatesToken_
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
edm::EDGetTokenT< std::vector< reco::Track > > recoTracksToken_
edm::EDGetTokenT< std::vector< TICLCandidate > > TICLCandidatesToken_
U second(std::pair< T, U > const &p)
edm::EDGetTokenT< ticl::RecoToSimCollectionSimTracksters > associatorMapRtSPUToken_
edm::EDGetTokenT< ticl::RecoToSimCollectionSimTracksters > associatorMapRtSToken_
Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge)
Definition: Common.h:44
edm::EDGetTokenT< ticl::SimToRecoCollectionSimTracksters > associatorMapStRToken_
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:141
edm::EDGetTokenT< std::vector< ticl::Trackster > > trackstersToken_
void bookCandidatesHistos(DQMStore::IBooker &ibook, Histograms &histograms, std::string baseDir) const
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
void fillCandidateHistos(const edm::Event &event, const Histograms &histograms, edm::Handle< ticl::TracksterCollection > simTrackstersCP_h) const
Definition: event.py:1