Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-16 05:06:46

0001 #include <numeric>
0002 #include <iomanip>
0003 #include <sstream>
0004 
0005 #include "Validation/HGCalValidation/interface/TICLCandidateValidator.h"
0006 #include "DataFormats/HGCalReco/interface/Common.h"
0007 
0008 TICLCandidateValidator::TICLCandidateValidator(edm::EDGetTokenT<std::vector<TICLCandidate>> ticlCandidates,
0009                                                edm::EDGetTokenT<std::vector<TICLCandidate>> simTICLCandidatesToken,
0010                                                edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken,
0011                                                edm::EDGetTokenT<std::vector<ticl::Trackster>> trackstersToken,
0012                                                edm::EDGetTokenT<ticl::TracksterToTracksterMap> associatorMapRtSToken,
0013                                                edm::EDGetTokenT<ticl::TracksterToTracksterMap> associatorMapStRToken,
0014                                                bool isTICLv5)
0015     : TICLCandidatesToken_(ticlCandidates),
0016       simTICLCandidatesToken_(simTICLCandidatesToken),
0017       recoTracksToken_(recoTracksToken),
0018       trackstersToken_(trackstersToken),
0019       associatorMapRtSToken_(associatorMapRtSToken),
0020       associatorMapStRToken_(associatorMapStRToken),
0021       isTICLv5_(isTICLv5) {}
0022 
0023 TICLCandidateValidator::~TICLCandidateValidator() {}
0024 
0025 void TICLCandidateValidator::bookCandidatesHistos(DQMStore::IBooker& ibook,
0026                                                   Histograms& histograms,
0027                                                   std::string baseDir) const {
0028   // book CAND histos
0029   histograms.h_tracksters_in_candidate =
0030       ibook.book1D("N of tracksters in candidate", "N of tracksters in candidate", 100, 0, 99);
0031   histograms.h_candidate_raw_energy =
0032       ibook.book1D("Candidates raw energy", "Candidates raw energy;E (GeV)", 250, 0, 250);
0033   histograms.h_candidate_regressed_energy =
0034       ibook.book1D("Candidates regressed energy", "Candidates regressed energy;E (GeV)", 250, 0, 250);
0035   histograms.h_candidate_pT = ibook.book1D("Candidates pT", "Candidates pT;p_{T}", 250, 0, 250);
0036   histograms.h_candidate_charge = ibook.book1D("Candidates charge", "Candidates charge;Charge", 3, -1.5, 1.5);
0037   histograms.h_candidate_pdgId = ibook.book1D("Candidates PDG Id", "Candidates PDG ID", 100, -220, 220);
0038   histograms.h_candidate_partType = ibook.book1D("Candidates type", "Candidates type", 9, -0.5, 8.5);
0039 
0040   // neutral: photon, pion, hadron
0041   const std::vector<std::string> neutrals{"photons", "neutral_pions", "neutral_hadrons"};
0042   for (long unsigned int i = 0; i < neutrals.size(); i++) {
0043     ibook.setCurrentFolder(baseDir + "/" + neutrals[i]);
0044 
0045     histograms.h_neut_tracksters_in_candidate.push_back(ibook.book1D("N of tracksters in candidate for " + neutrals[i],
0046                                                                      "N of tracksters in candidate for " + neutrals[i],
0047                                                                      100,
0048                                                                      0,
0049                                                                      99));
0050     histograms.h_neut_candidate_regressed_energy.push_back(ibook.book1D(
0051         neutrals[i] + "candidates regressed energy", neutrals[i] + " candidates regressed energy;E (GeV)", 250, 0, 250));
0052     histograms.h_neut_candidate_charge.push_back(
0053         ibook.book1D(neutrals[i] + " candidates charge", neutrals[i] + " candidates charge;Charge", 3, -1.5, 1.5));
0054     histograms.h_neut_candidate_pdgId.push_back(
0055         ibook.book1D(neutrals[i] + " candidates PDG Id", neutrals[i] + " candidates PDG ID", 100, -220, 220));
0056     histograms.h_neut_candidate_partType.push_back(
0057         ibook.book1D(neutrals[i] + " candidates type", neutrals[i] + " candidates type", 9, -0.5, 8.5));
0058 
0059     histograms.h_den_fake_neut_energy_candidate.push_back(
0060         ibook.book1D("den_fake_cand_vs_energy_" + neutrals[i], neutrals[i] + " candidates energy;E (GeV)", 50, 0, 250));
0061     histograms.h_num_fake_neut_energy_candidate_pdgId.push_back(ibook.book1D(
0062         "num_fake_pid_cand_vs_energy_" + neutrals[i], neutrals[i] + " PID fake vs energy;E (GeV)", 50, 0, 250));
0063     histograms.h_num_fake_neut_energy_candidate_energy.push_back(
0064         ibook.book1D("num_fake_energy_cand_vs_energy_" + neutrals[i],
0065                      neutrals[i] + " PID and energy fake vs energy;E (GeV)",
0066                      50,
0067                      0,
0068                      250));
0069     histograms.h_den_fake_neut_pt_candidate.push_back(
0070         ibook.book1D("den_fake_cand_vs_pt_" + neutrals[i], neutrals[i] + " candidates pT;p_{T} (GeV)", 50, 0, 250));
0071     histograms.h_num_fake_neut_pt_candidate_pdgId.push_back(ibook.book1D(
0072         "num_fake_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " PID fake vs pT;p_{T} (GeV)", 50, 0, 250));
0073     histograms.h_num_fake_neut_pt_candidate_energy.push_back(
0074         ibook.book1D("num_fake_energy_cand_vs_pt_" + neutrals[i],
0075                      neutrals[i] + " PID and energy fake vs pT;p_{T} (GeV)",
0076                      50,
0077                      0,
0078                      250));
0079     histograms.h_den_fake_neut_eta_candidate.push_back(
0080         ibook.book1D("den_fake_cand_vs_eta_" + neutrals[i], neutrals[i] + " candidates eta;#eta (GeV)", 50, -3, 3));
0081     histograms.h_num_fake_neut_eta_candidate_pdgId.push_back(ibook.book1D(
0082         "num_fake_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " PID fake vs eta;#eta (GeV)", 50, -3, 3));
0083     histograms.h_num_fake_neut_eta_candidate_energy.push_back(
0084         ibook.book1D("num_fake_energy_cand_vs_eta_" + neutrals[i],
0085                      neutrals[i] + " PID and energy fake vs eta;#eta (GeV)",
0086                      50,
0087                      -3,
0088                      3));
0089     histograms.h_den_fake_neut_phi_candidate.push_back(ibook.book1D(
0090         "den_fake_cand_vs_phi_" + neutrals[i], neutrals[i] + " candidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
0091     histograms.h_num_fake_neut_phi_candidate_pdgId.push_back(ibook.book1D(
0092         "num_fake_pid_cand_vs_phi_" + neutrals[i], neutrals[i] + " PID fake vs phi;#phi (GeV)", 50, -3.14159, 3.14159));
0093     histograms.h_num_fake_neut_phi_candidate_energy.push_back(
0094         ibook.book1D("num_fake_energy_cand_vs_phi_" + neutrals[i],
0095                      neutrals[i] + " PID and energy fake vs phi;#phi (GeV)",
0096                      50,
0097                      -3.14159,
0098                      3.14159));
0099 
0100     histograms.h_den_neut_energy_candidate.push_back(
0101         ibook.book1D("den_cand_vs_energy_" + neutrals[i], neutrals[i] + " simCandidates energy;E (GeV)", 50, 0, 250));
0102     histograms.h_num_neut_energy_candidate_pdgId.push_back(
0103         ibook.book1D("num_pid_cand_vs_energy_" + neutrals[i],
0104                      neutrals[i] + " track and PID efficiency vs energy;E (GeV)",
0105                      50,
0106                      0,
0107                      250));
0108     histograms.h_num_neut_energy_candidate_energy.push_back(
0109         ibook.book1D("num_energy_cand_vs_energy_" + neutrals[i],
0110                      neutrals[i] + " track, PID and energy efficiency vs energy;E (GeV)",
0111                      50,
0112                      0,
0113                      250));
0114     histograms.h_den_neut_pt_candidate.push_back(
0115         ibook.book1D("den_cand_vs_pt_" + neutrals[i], neutrals[i] + " simCandidates pT;p_{T} (GeV)", 50, 0, 250));
0116     histograms.h_num_neut_pt_candidate_pdgId.push_back(ibook.book1D(
0117         "num_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
0118     histograms.h_num_neut_pt_candidate_energy.push_back(
0119         ibook.book1D("num_energy_cand_vs_pt_" + neutrals[i],
0120                      neutrals[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)",
0121                      50,
0122                      0,
0123                      250));
0124     histograms.h_den_neut_eta_candidate.push_back(
0125         ibook.book1D("den_cand_vs_eta_" + neutrals[i], neutrals[i] + " simCandidates eta;#eta (GeV)", 50, -3, 3));
0126     histograms.h_num_neut_eta_candidate_pdgId.push_back(ibook.book1D(
0127         "num_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " track and PID efficiency vs eta;#eta (GeV)", 50, -3, 3));
0128     histograms.h_num_neut_eta_candidate_energy.push_back(
0129         ibook.book1D("num_energy_cand_vs_eta_" + neutrals[i],
0130                      neutrals[i] + " track, PID and energy efficiency vs eta;#eta (GeV)",
0131                      50,
0132                      -3,
0133                      3));
0134     histograms.h_den_neut_phi_candidate.push_back(ibook.book1D(
0135         "den_cand_vs_phi_" + neutrals[i], neutrals[i] + " simCandidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
0136     histograms.h_num_neut_phi_candidate_pdgId.push_back(
0137         ibook.book1D("num_pid_cand_vs_phi_" + neutrals[i],
0138                      neutrals[i] + " track and PID efficiency vs phi;#phi (GeV)",
0139                      50,
0140                      -3.14159,
0141                      3.14159));
0142     histograms.h_num_neut_phi_candidate_energy.push_back(
0143         ibook.book1D("num_energy_cand_vs_phi_" + neutrals[i],
0144                      neutrals[i] + " track, PID and energy efficiency vs phi;#phi (GeV)",
0145                      50,
0146                      -3.14159,
0147                      3.14159));
0148   }
0149   // charged: electron, muon, hadron
0150   const std::vector<std::string> charged{"electrons", "muons", "charged_hadrons"};
0151   for (long unsigned int i = 0; i < charged.size(); i++) {
0152     ibook.setCurrentFolder(baseDir + "/" + charged[i]);
0153 
0154     histograms.h_chg_tracksters_in_candidate.push_back(ibook.book1D(
0155         "N of tracksters in candidate for " + charged[i], "N of tracksters in candidate for " + charged[i], 100, 0, 99));
0156     histograms.h_chg_candidate_regressed_energy.push_back(ibook.book1D(
0157         charged[i] + "candidates regressed energy", charged[i] + " candidates regressed energy;E (GeV)", 250, 0, 250));
0158     histograms.h_chg_candidate_charge.push_back(
0159         ibook.book1D(charged[i] + " candidates charge", charged[i] + " candidates charge;Charge", 3, -1.5, 1.5));
0160     histograms.h_chg_candidate_pdgId.push_back(
0161         ibook.book1D(charged[i] + " candidates PDG Id", charged[i] + " candidates PDG ID", 100, -220, 220));
0162     histograms.h_chg_candidate_partType.push_back(
0163         ibook.book1D(charged[i] + " candidates type", charged[i] + " candidates type", 9, -0.5, 8.5));
0164 
0165     histograms.h_den_fake_chg_energy_candidate.push_back(
0166         ibook.book1D("den_fake_cand_vs_energy_" + charged[i], charged[i] + " candidates energy;E (GeV)", 50, 0, 250));
0167     histograms.h_num_fake_chg_energy_candidate_track.push_back(ibook.book1D(
0168         "num_fake_track_cand_vs_energy_" + charged[i], charged[i] + " track fake vs energy;E (GeV)", 50, 0, 250));
0169     histograms.h_num_fake_chg_energy_candidate_pdgId.push_back(ibook.book1D(
0170         "num_fake_pid_cand_vs_energy_" + charged[i], charged[i] + " track and PID fake vs energy;E (GeV)", 50, 0, 250));
0171     histograms.h_num_fake_chg_energy_candidate_energy.push_back(
0172         ibook.book1D("num_fake_energy_cand_vs_energy_" + charged[i],
0173                      charged[i] + " track, PID and energy fake vs energy;E (GeV)",
0174                      50,
0175                      0,
0176                      250));
0177     histograms.h_den_fake_chg_pt_candidate.push_back(
0178         ibook.book1D("den_fake_cand_vs_pt_" + charged[i], charged[i] + " candidates pT;p_{T} (GeV)", 50, 0, 250));
0179     histograms.h_num_fake_chg_pt_candidate_track.push_back(ibook.book1D(
0180         "num_fake_track_cand_vs_pt_" + charged[i], charged[i] + " track fake vs pT;p_{T} (GeV)", 50, 0, 250));
0181     histograms.h_num_fake_chg_pt_candidate_pdgId.push_back(ibook.book1D(
0182         "num_fake_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID fake vs pT;p_{T} (GeV)", 50, 0, 250));
0183     histograms.h_num_fake_chg_pt_candidate_energy.push_back(
0184         ibook.book1D("num_fake_energy_cand_vs_pt_" + charged[i],
0185                      charged[i] + " track, PID and energy fake vs pT;p_{T} (GeV)",
0186                      50,
0187                      0,
0188                      250));
0189     histograms.h_den_fake_chg_eta_candidate.push_back(
0190         ibook.book1D("den_fake_cand_vs_eta_" + charged[i], charged[i] + " candidates eta;#eta (GeV)", 50, -3, 3));
0191     histograms.h_num_fake_chg_eta_candidate_track.push_back(ibook.book1D(
0192         "num_fake_track_cand_vs_eta_" + charged[i], charged[i] + " track fake vs eta;#eta (GeV)", 50, -3, 3));
0193     histograms.h_num_fake_chg_eta_candidate_pdgId.push_back(ibook.book1D(
0194         "num_fake_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID fake vs eta;#eta (GeV)", 50, -3, 3));
0195     histograms.h_num_fake_chg_eta_candidate_energy.push_back(
0196         ibook.book1D("num_fake_energy_cand_vs_eta_" + charged[i],
0197                      charged[i] + " track, PID and energy fake vs eta;#eta (GeV)",
0198                      50,
0199                      -3,
0200                      3));
0201     histograms.h_den_fake_chg_phi_candidate.push_back(ibook.book1D(
0202         "den_fake_cand_vs_phi_" + charged[i], charged[i] + " candidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
0203     histograms.h_num_fake_chg_phi_candidate_track.push_back(ibook.book1D("num_fake_track_cand_vs_phi_" + charged[i],
0204                                                                          charged[i] + " track fake vs phi;#phi (GeV)",
0205                                                                          50,
0206                                                                          -3.14159,
0207                                                                          3.14159));
0208     histograms.h_num_fake_chg_phi_candidate_pdgId.push_back(
0209         ibook.book1D("num_fake_pid_cand_vs_phi_" + charged[i],
0210                      charged[i] + " track and PID fake vs phi;#phi (GeV)",
0211                      50,
0212                      -3.14159,
0213                      3.14159));
0214     histograms.h_num_fake_chg_phi_candidate_energy.push_back(
0215         ibook.book1D("num_fake_energy_cand_vs_phi_" + charged[i],
0216                      charged[i] + " track, PID and energy fake vs phi;#phi (GeV)",
0217                      50,
0218                      -3.14159,
0219                      3.14159));
0220 
0221     histograms.h_den_chg_energy_candidate.push_back(
0222         ibook.book1D("den_cand_vs_energy_" + charged[i], charged[i] + " simCandidates energy;E (GeV)", 50, 0, 250));
0223     histograms.h_num_chg_energy_candidate_track.push_back(ibook.book1D(
0224         "num_track_cand_vs_energy_" + charged[i], charged[i] + " track efficiency vs energy;E (GeV)", 50, 0, 250));
0225     histograms.h_num_chg_energy_candidate_pdgId.push_back(ibook.book1D(
0226         "num_pid_cand_vs_energy_" + charged[i], charged[i] + " track and PID efficiency vs energy;E (GeV)", 50, 0, 250));
0227     histograms.h_num_chg_energy_candidate_energy.push_back(
0228         ibook.book1D("num_energy_cand_vs_energy_" + charged[i],
0229                      charged[i] + " track, PID and energy efficiency vs energy;E (GeV)",
0230                      50,
0231                      0,
0232                      250));
0233     histograms.h_den_chg_pt_candidate.push_back(
0234         ibook.book1D("den_cand_vs_pt_" + charged[i], charged[i] + " simCandidates pT;p_{T} (GeV)", 50, 0, 250));
0235     histograms.h_num_chg_pt_candidate_track.push_back(ibook.book1D(
0236         "num_track_cand_vs_pt_" + charged[i], charged[i] + " track efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
0237     histograms.h_num_chg_pt_candidate_pdgId.push_back(ibook.book1D(
0238         "num_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 50, 0, 250));
0239     histograms.h_num_chg_pt_candidate_energy.push_back(
0240         ibook.book1D("num_energy_cand_vs_pt_" + charged[i],
0241                      charged[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)",
0242                      50,
0243                      0,
0244                      250));
0245     histograms.h_den_chg_eta_candidate.push_back(
0246         ibook.book1D("den_cand_vs_eta_" + charged[i], charged[i] + " simCandidates eta;#eta (GeV)", 50, -3, 3));
0247     histograms.h_num_chg_eta_candidate_track.push_back(ibook.book1D(
0248         "num_track_cand_vs_eta_" + charged[i], charged[i] + " track efficiency vs eta;#eta (GeV)", 50, -3, 3));
0249     histograms.h_num_chg_eta_candidate_pdgId.push_back(ibook.book1D(
0250         "num_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID efficiency vs eta;#eta (GeV)", 50, -3, 3));
0251     histograms.h_num_chg_eta_candidate_energy.push_back(
0252         ibook.book1D("num_energy_cand_vs_eta_" + charged[i],
0253                      charged[i] + " track, PID and energy efficiency vs eta;#eta (GeV)",
0254                      50,
0255                      -3,
0256                      3));
0257     histograms.h_den_chg_phi_candidate.push_back(ibook.book1D(
0258         "den_cand_vs_phi_" + charged[i], charged[i] + " simCandidates phi;#phi (GeV)", 50, -3.14159, 3.14159));
0259     histograms.h_num_chg_phi_candidate_track.push_back(ibook.book1D("num_track_cand_vs_phi_" + charged[i],
0260                                                                     charged[i] + " track efficiency vs phi;#phi (GeV)",
0261                                                                     50,
0262                                                                     -3.14159,
0263                                                                     3.14159));
0264     histograms.h_num_chg_phi_candidate_pdgId.push_back(
0265         ibook.book1D("num_pid_cand_vs_phi_" + charged[i],
0266                      charged[i] + " track and PID efficiency vs phi;#phi (GeV)",
0267                      50,
0268                      -3.14159,
0269                      3.14159));
0270     histograms.h_num_chg_phi_candidate_energy.push_back(
0271         ibook.book1D("num_energy_cand_vs_phi_" + charged[i],
0272                      charged[i] + " track, PID and energy efficiency vs phi;#phi (GeV)",
0273                      50,
0274                      -3.14159,
0275                      3.14159));
0276   }
0277 }
0278 
0279 void TICLCandidateValidator::fillCandidateHistos(const edm::Event& event,
0280                                                  const Histograms& histograms,
0281                                                  edm::Handle<ticl::TracksterCollection> simTrackstersCP_h) const {
0282   auto TICLCandidates = event.get(TICLCandidatesToken_);
0283 
0284   edm::Handle<std::vector<TICLCandidate>> simTICLCandidates_h;
0285   event.getByToken(simTICLCandidatesToken_, simTICLCandidates_h);
0286   auto simTICLCandidates = *simTICLCandidates_h;
0287 
0288   edm::Handle<std::vector<reco::Track>> recoTracks_h;
0289   event.getByToken(recoTracksToken_, recoTracks_h);
0290   auto recoTracks = *recoTracks_h;
0291 
0292   edm::Handle<std::vector<ticl::Trackster>> Tracksters_h;
0293   event.getByToken(trackstersToken_, Tracksters_h);
0294   auto trackstersMerged = *Tracksters_h;
0295 
0296   edm::Handle<ticl::TracksterToTracksterMap> mergeTsRecoToSim_h;
0297   event.getByToken(associatorMapRtSToken_, mergeTsRecoToSim_h);
0298   auto const& mergeTsRecoToSimMap = *mergeTsRecoToSim_h;
0299 
0300   edm::Handle<ticl::TracksterToTracksterMap> mergeTsSimToReco_h;
0301   event.getByToken(associatorMapStRToken_, mergeTsSimToReco_h);
0302   auto const& mergeTsSimToRecoMap = *mergeTsSimToReco_h;
0303 
0304   // candidates plots
0305   for (const auto& cand : TICLCandidates) {
0306     histograms.h_tracksters_in_candidate->Fill(cand.tracksters().size());
0307     histograms.h_candidate_raw_energy->Fill(cand.rawEnergy());
0308     histograms.h_candidate_regressed_energy->Fill(cand.energy());
0309     histograms.h_candidate_pT->Fill(cand.pt());
0310     histograms.h_candidate_charge->Fill(cand.charge());
0311     histograms.h_candidate_pdgId->Fill(cand.pdgId());
0312     const auto& arr = cand.idProbabilities();
0313     histograms.h_candidate_partType->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
0314   }
0315 
0316   std::vector<int> chargedCandidates;
0317   std::vector<int> neutralCandidates;
0318   chargedCandidates.reserve(simTICLCandidates.size());
0319   neutralCandidates.reserve(simTICLCandidates.size());
0320 
0321   for (size_t i = 0; i < simTICLCandidates.size(); ++i) {
0322     const auto& simCand = simTICLCandidates[i];
0323     const auto particleType = ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), simCand.charge());
0324     if (particleType == ticl::Trackster::ParticleType::electron or
0325         particleType == ticl::Trackster::ParticleType::muon or
0326         particleType == ticl::Trackster::ParticleType::charged_hadron)
0327       chargedCandidates.emplace_back(i);
0328     else if (particleType == ticl::Trackster::ParticleType::photon or
0329              particleType == ticl::Trackster::ParticleType::neutral_pion or
0330              particleType == ticl::Trackster::ParticleType::neutral_hadron)
0331       neutralCandidates.emplace_back(i);
0332     // should consider also unknown ?
0333   }
0334 
0335   chargedCandidates.shrink_to_fit();
0336   neutralCandidates.shrink_to_fit();
0337 
0338   for (const auto i : chargedCandidates) {
0339     const auto& simCand = simTICLCandidates[i];
0340     auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 1)));
0341     /* 11 (type 1) becomes 0
0342      * 13 (type 2) becomes 1
0343      * 211 (type 4) becomes 2
0344      */
0345     int32_t simCandTrackIdx = -1;
0346     if (simCand.trackPtr().get() != nullptr)
0347       simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
0348     else {
0349       // no reco track, but simCand is charged
0350       continue;
0351     }
0352     if (simCand.trackPtr().get()->pt() < 1 or simCand.trackPtr().get()->missingOuterHits() > 5 or
0353         not simCand.trackPtr().get()->quality(reco::TrackBase::highPurity))
0354       continue;
0355 
0356     // +1 to all denominators
0357     histograms.h_den_chg_energy_candidate[index]->Fill(simCand.rawEnergy());
0358     histograms.h_den_chg_pt_candidate[index]->Fill(simCand.pt());
0359     histograms.h_den_chg_eta_candidate[index]->Fill(simCand.eta());
0360     histograms.h_den_chg_phi_candidate[index]->Fill(simCand.phi());
0361 
0362     int32_t cand_idx = -1;
0363     float shared_energy = 0.;
0364     const auto ts_vec = mergeTsSimToRecoMap[i];
0365     if (!ts_vec.empty()) {
0366       auto min_elem =
0367           std::min_element(ts_vec.begin(), ts_vec.end(), [](auto const& ts1_id_pair, auto const& ts2_id_pair) {
0368             return ts1_id_pair.second.second < ts2_id_pair.second.second;
0369           });
0370       shared_energy = min_elem->second.first;
0371       cand_idx = min_elem->first;
0372     }
0373     // no reco associated to sim
0374     if (cand_idx == -1)
0375       continue;
0376 
0377     auto& recoCand = TICLCandidates[cand_idx];
0378     if (isTICLv5_) {
0379       // cand_idx is the tsMerge index, find the ts in the candidates collection
0380       auto const tsPtr = edm::Ptr<ticl::Trackster>(Tracksters_h, cand_idx);
0381       auto cand_it = std::find_if(TICLCandidates.begin(), TICLCandidates.end(), [tsPtr](TICLCandidate const& cand) {
0382         if (!cand.tracksters().empty())
0383           return cand.tracksters()[0] == tsPtr;
0384         else
0385           return false;
0386       });
0387       if (cand_it != TICLCandidates.end())
0388         recoCand = *cand_it;
0389       else
0390         continue;
0391     }
0392 
0393     if (recoCand.trackPtr().get() != nullptr) {
0394       const auto candTrackIdx = recoCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
0395       if (simCandTrackIdx == candTrackIdx) {
0396         // +1 to track num
0397         histograms.h_num_chg_energy_candidate_track[index]->Fill(simCand.rawEnergy());
0398         histograms.h_num_chg_pt_candidate_track[index]->Fill(simCand.pt());
0399         histograms.h_num_chg_eta_candidate_track[index]->Fill(simCand.eta());
0400         histograms.h_num_chg_phi_candidate_track[index]->Fill(simCand.phi());
0401       } else {
0402         continue;
0403       }
0404     } else {
0405       continue;
0406     }
0407 
0408     //step 2: PID
0409     if (simCand.pdgId() == recoCand.pdgId()) {
0410       // +1 to num pdg id
0411       histograms.h_num_chg_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy());
0412       histograms.h_num_chg_pt_candidate_pdgId[index]->Fill(simCand.pt());
0413       histograms.h_num_chg_eta_candidate_pdgId[index]->Fill(simCand.eta());
0414       histograms.h_num_chg_phi_candidate_pdgId[index]->Fill(simCand.phi());
0415 
0416       //step 3: energy
0417       if (shared_energy / simCand.rawEnergy() > 0.5) {
0418         // +1 to ene num
0419         histograms.h_num_chg_energy_candidate_energy[index]->Fill(simCand.rawEnergy());
0420         histograms.h_num_chg_pt_candidate_energy[index]->Fill(simCand.pt());
0421         histograms.h_num_chg_eta_candidate_energy[index]->Fill(simCand.eta());
0422         histograms.h_num_chg_phi_candidate_energy[index]->Fill(simCand.phi());
0423       }
0424     }
0425   }
0426 
0427   for (const auto i : neutralCandidates) {
0428     const auto& simCand = simTICLCandidates[i];
0429     auto index = int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 0)) / 2;
0430     /* 22 (type 0) becomes 0
0431      * 111 (type 3) becomes 1
0432      * 130 (type 5) becomes 2
0433      */
0434     histograms.h_den_neut_energy_candidate[index]->Fill(simCand.rawEnergy());
0435     histograms.h_den_neut_pt_candidate[index]->Fill(simCand.pt());
0436     histograms.h_den_neut_eta_candidate[index]->Fill(simCand.eta());
0437     histograms.h_den_neut_phi_candidate[index]->Fill(simCand.phi());
0438 
0439     int32_t cand_idx = -1;
0440     float shared_energy = 0.;
0441     const auto ts_vec = mergeTsSimToRecoMap[i];
0442     if (!ts_vec.empty()) {
0443       auto min_elem =
0444           std::min_element(ts_vec.begin(), ts_vec.end(), [](auto const& ts1_id_pair, auto const& ts2_id_pair) {
0445             return ts1_id_pair.second.second < ts2_id_pair.second.second;
0446           });
0447       shared_energy = min_elem->second.first;
0448       cand_idx = min_elem->first;
0449     }
0450 
0451     // no reco associated to sim
0452     if (cand_idx == -1)
0453       continue;
0454 
0455     auto& recoCand = TICLCandidates[cand_idx];
0456     if (isTICLv5_) {
0457       // cand_idx is the tsMerge index, find the ts in the candidates collection
0458       auto const tsPtr = edm::Ptr<ticl::Trackster>(Tracksters_h, cand_idx);
0459       auto cand_it = std::find_if(TICLCandidates.begin(), TICLCandidates.end(), [tsPtr](TICLCandidate const& cand) {
0460         if (!cand.tracksters().empty())
0461           return cand.tracksters()[0] == tsPtr;
0462         else
0463           return false;
0464       });
0465       if (cand_it != TICLCandidates.end())
0466         recoCand = *cand_it;
0467       else
0468         continue;
0469     }
0470 
0471     if (recoCand.trackPtr().get() != nullptr)
0472       continue;
0473 
0474     //step 2: PID
0475     if (simCand.pdgId() == recoCand.pdgId()) {
0476       // +1 to num pdg id
0477       histograms.h_num_neut_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy());
0478       histograms.h_num_neut_pt_candidate_pdgId[index]->Fill(simCand.pt());
0479       histograms.h_num_neut_eta_candidate_pdgId[index]->Fill(simCand.eta());
0480       histograms.h_num_neut_phi_candidate_pdgId[index]->Fill(simCand.phi());
0481 
0482       //step 3: energy
0483       if (shared_energy / simCand.rawEnergy() > 0.5) {
0484         // +1 to ene num
0485         histograms.h_num_neut_energy_candidate_energy[index]->Fill(simCand.rawEnergy());
0486         histograms.h_num_neut_pt_candidate_energy[index]->Fill(simCand.pt());
0487         histograms.h_num_neut_eta_candidate_energy[index]->Fill(simCand.eta());
0488         histograms.h_num_neut_phi_candidate_energy[index]->Fill(simCand.phi());
0489       }
0490     }
0491   }
0492 
0493   // FAKE rate
0494   chargedCandidates.clear();
0495   neutralCandidates.clear();
0496   chargedCandidates.reserve(TICLCandidates.size());
0497   neutralCandidates.reserve(TICLCandidates.size());
0498 
0499   auto isCharged = [](int pdgId) {
0500     pdgId = std::abs(pdgId);
0501     return (pdgId == 11 or pdgId == 211 or pdgId == 13);
0502   };
0503 
0504   for (size_t i = 0; i < TICLCandidates.size(); ++i) {
0505     const auto& cand = TICLCandidates[i];
0506     const auto& charged = isCharged(cand.pdgId());
0507     if (charged)
0508       chargedCandidates.emplace_back(i);
0509     else
0510       neutralCandidates.emplace_back(i);
0511 
0512     // should consider also unknown ?
0513   }
0514 
0515   chargedCandidates.shrink_to_fit();
0516   neutralCandidates.shrink_to_fit();
0517 
0518   // loop on charged
0519   for (const auto i : chargedCandidates) {
0520     const auto& cand = TICLCandidates[i];
0521     auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 1)));
0522     /* 11 (type 1) becomes 0
0523      * 13 (type 2) becomes 1
0524      * 211 (type 4) becomes 2
0525      */
0526     int32_t candTrackIdx = -1;
0527     candTrackIdx = cand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
0528 
0529     if (cand.tracksters().empty())
0530       continue;
0531 
0532     // i is the candidate idx == ts idx only in v4, find ts_idx in v5
0533     auto mergeTs_id = i;
0534     if (isTICLv5_) {
0535       mergeTs_id = cand.tracksters()[0].get() - edm::Ptr<ticl::Trackster>(Tracksters_h, 0).get();
0536     }
0537 
0538     // +1 to all denominators
0539     histograms.h_den_fake_chg_energy_candidate[index]->Fill(cand.rawEnergy());
0540     histograms.h_den_fake_chg_pt_candidate[index]->Fill(cand.pt());
0541     histograms.h_den_fake_chg_eta_candidate[index]->Fill(cand.eta());
0542     histograms.h_den_fake_chg_phi_candidate[index]->Fill(cand.phi());
0543 
0544     histograms.h_chg_tracksters_in_candidate[index]->Fill(cand.tracksters().size());
0545     histograms.h_chg_candidate_regressed_energy[index]->Fill(cand.energy());
0546     histograms.h_chg_candidate_charge[index]->Fill(cand.charge());
0547     histograms.h_chg_candidate_pdgId[index]->Fill(cand.pdgId());
0548     const auto& arr = cand.idProbabilities();
0549     histograms.h_chg_candidate_partType[index]->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
0550 
0551     int32_t simCand_idx = -1;
0552     const auto sts_vec = mergeTsRecoToSimMap[mergeTs_id];
0553     float shared_energy = 0.;
0554     // search for reco cand associated
0555     if (!sts_vec.empty()) {
0556       auto min_elem =
0557           std::min_element(sts_vec.begin(), sts_vec.end(), [](auto const& sts1_id_pair, auto const& sts2_id_pair) {
0558             return sts1_id_pair.second.second < sts2_id_pair.second.second;
0559           });
0560       shared_energy = min_elem->second.first;
0561       simCand_idx = min_elem->first;
0562     }
0563 
0564     if (simCand_idx == -1)
0565       continue;
0566 
0567     const auto& simCand = simTICLCandidates[simCand_idx];
0568     if (simCand.trackPtr().get() != nullptr) {
0569       const auto simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr<reco::Track>(recoTracks_h, 0).get();
0570       if (simCandTrackIdx != candTrackIdx) {
0571         // fake += 1
0572         histograms.h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy());
0573         histograms.h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt());
0574         histograms.h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta());
0575         histograms.h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi());
0576         continue;
0577       }
0578     } else {
0579       // fake += 1
0580       histograms.h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy());
0581       histograms.h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt());
0582       histograms.h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta());
0583       histograms.h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi());
0584       continue;
0585     }
0586 
0587     //step 2: PID
0588     if (simCand.pdgId() != cand.pdgId()) {
0589       // +1 to num fake pdg id
0590       histograms.h_num_fake_chg_energy_candidate_pdgId[index]->Fill(cand.rawEnergy());
0591       histograms.h_num_fake_chg_pt_candidate_pdgId[index]->Fill(cand.pt());
0592       histograms.h_num_fake_chg_eta_candidate_pdgId[index]->Fill(cand.eta());
0593       histograms.h_num_fake_chg_phi_candidate_pdgId[index]->Fill(cand.phi());
0594       continue;
0595     }
0596 
0597     //step 3: energy
0598     if (shared_energy / simCand.rawEnergy() < 0.5) {
0599       // +1 to ene num
0600       histograms.h_num_fake_chg_energy_candidate_energy[index]->Fill(cand.rawEnergy());
0601       histograms.h_num_fake_chg_pt_candidate_energy[index]->Fill(cand.pt());
0602       histograms.h_num_fake_chg_eta_candidate_energy[index]->Fill(cand.eta());
0603       histograms.h_num_fake_chg_phi_candidate_energy[index]->Fill(cand.phi());
0604     }
0605   }
0606   // loop on neutrals
0607   for (const auto i : neutralCandidates) {
0608     const auto& cand = TICLCandidates[i];
0609     auto index = int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 0)) / 2;
0610     /* 22 (type 0) becomes 0
0611      * 111 (type 3) becomes 1
0612      * 130 (type 5) becomes 2
0613      */
0614 
0615     if (cand.tracksters().empty())
0616       continue;
0617 
0618     // i is the candidate idx == ts idx only in v4, find ts_idx in v5
0619     auto mergeTs_id = i;
0620     if (isTICLv5_) {
0621       mergeTs_id = cand.tracksters()[0].get() - edm::Ptr<ticl::Trackster>(Tracksters_h, 0).get();
0622     }
0623 
0624     // +1 to all denominators
0625     histograms.h_den_fake_neut_energy_candidate[index]->Fill(cand.rawEnergy());
0626     histograms.h_den_fake_neut_pt_candidate[index]->Fill(cand.pt());
0627     histograms.h_den_fake_neut_eta_candidate[index]->Fill(cand.eta());
0628     histograms.h_den_fake_neut_phi_candidate[index]->Fill(cand.phi());
0629 
0630     histograms.h_neut_tracksters_in_candidate[index]->Fill(cand.tracksters().size());
0631     histograms.h_neut_candidate_regressed_energy[index]->Fill(cand.energy());
0632     histograms.h_neut_candidate_charge[index]->Fill(cand.charge());
0633     histograms.h_neut_candidate_pdgId[index]->Fill(cand.pdgId());
0634     const auto& arr = cand.idProbabilities();
0635     histograms.h_neut_candidate_partType[index]->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin());
0636 
0637     int32_t simCand_idx = -1;
0638     const auto sts_vec = mergeTsRecoToSimMap[mergeTs_id];
0639     float shared_energy = 0.;
0640     // search for reco cand associated
0641     if (!sts_vec.empty()) {
0642       auto min_elem =
0643           std::min_element(sts_vec.begin(), sts_vec.end(), [](auto const& sts1_id_pair, auto const& sts2_id_pair) {
0644             return sts1_id_pair.second.second < sts2_id_pair.second.second;
0645           });
0646       shared_energy = min_elem->second.first;
0647       simCand_idx = min_elem->first;
0648     }
0649 
0650     if (simCand_idx == -1)
0651       continue;
0652 
0653     const auto& simCand = simTICLCandidates[simCand_idx];
0654 
0655     //step 2: PID
0656     if (simCand.pdgId() != cand.pdgId()) {
0657       // +1 to num fake pdg id
0658       histograms.h_num_fake_neut_energy_candidate_pdgId[index]->Fill(cand.rawEnergy());
0659       histograms.h_num_fake_neut_pt_candidate_pdgId[index]->Fill(cand.pt());
0660       histograms.h_num_fake_neut_eta_candidate_pdgId[index]->Fill(cand.eta());
0661       histograms.h_num_fake_neut_phi_candidate_pdgId[index]->Fill(cand.phi());
0662       continue;
0663     }
0664 
0665     //step 3: energy
0666     if (shared_energy / simCand.rawEnergy() < 0.5) {
0667       // +1 to ene num
0668       histograms.h_num_fake_neut_energy_candidate_energy[index]->Fill(cand.rawEnergy());
0669       histograms.h_num_fake_neut_pt_candidate_energy[index]->Fill(cand.pt());
0670       histograms.h_num_fake_neut_eta_candidate_energy[index]->Fill(cand.eta());
0671       histograms.h_num_fake_neut_phi_candidate_energy[index]->Fill(cand.phi());
0672     }
0673   }
0674 }