From a2b82587bd50e560e2cab66b22003cd5db4a90d0 Mon Sep 17 00:00:00 2001 From: Aurora Perego <84700885+AuroraPerego@users.noreply.github.com> Date: Thu, 21 Dec 2023 13:27:54 +0100 Subject: [PATCH] TICLCandidate validation (#21) * changes to simTICLCandidates pdgId - keep original pdgId for charged simTICLCandidates w/o reco track - add check for pi0 pdgId * add ticlCandidates validation --- .../TICL/plugins/SimTrackstersProducer.cc | 13 +- .../interface/HGCalValidator.h | 4 + .../interface/TICLCandidateValidator.h | 118 ++++ .../HGCalValidation/plugins/HGCalValidator.cc | 37 + .../python/HGCalPostProcessor_cff.py | 5 +- .../python/HGCalValidator_cfi.py | 10 + .../python/PostProcessorHGCAL_cfi.py | 36 + .../src/TICLCandidateValidator.cc | 655 ++++++++++++++++++ 8 files changed, 872 insertions(+), 6 deletions(-) create mode 100644 Validation/HGCalValidation/interface/TICLCandidateValidator.h create mode 100644 Validation/HGCalValidation/src/TICLCandidateValidator.cc diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index c062f36a48fa8..1bd563923a0ce 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -459,9 +459,7 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) auto pdgId = cp.pdgId(); auto charge = cp.charge(); - if (cand.trackPtr().isNonnull() and charge == 0) { - } - if (cand.trackPtr().isNonnull() and charge != 0) { + if (cand.trackPtr().isNonnull()) { auto const& track = cand.trackPtr().get(); if (std::abs(pdgId) == 13) { cand.setPdgId(pdgId); @@ -475,7 +473,14 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) regressedEnergy); cand.setP4(p4); } else { // neutral candidates - cand.setPdgId(isHad(pdgId) ? 130 : 22); + // a neutral candidate with a charged CaloParticle is charged without a reco track associated with it + // set the charge = 0, but keep the real pdgId to keep track of that + if (charge != 0) + cand.setPdgId(isHad(pdgId) ? 211 : 11); + else if (pdgId == 111) + cand.setPdgId(pdgId); + else + cand.setPdgId(isHad(pdgId) ? 130 : 22); cand.setCharge(0); auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1); diff --git a/Validation/HGCalValidation/interface/HGCalValidator.h b/Validation/HGCalValidation/interface/HGCalValidator.h index 48fc3dc52f7b0..f2c8901c4bca9 100644 --- a/Validation/HGCalValidation/interface/HGCalValidator.h +++ b/Validation/HGCalValidation/interface/HGCalValidator.h @@ -22,6 +22,7 @@ #include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h" +#include "Validation/HGCalValidation/interface/TICLCandidateValidator.h" #include "Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h" #include "Validation/HGCalValidation/interface/CaloParticleSelector.h" #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h" @@ -75,6 +76,8 @@ class HGCalValidator : public DQMGlobalEDAnalyzer { const bool doTrackstersPlots_; std::string label_TS_, label_TSToCPLinking_, label_TSToSTSPR_; std::vector label_clustersmask; + const bool doCandidatesPlots_; + edm::InputTag label_candidates_; const edm::FileInPath cummatbudinxo_; std::vector> labelToken; @@ -96,6 +99,7 @@ class HGCalValidator : public DQMGlobalEDAnalyzer { std::unique_ptr histoProducerAlgo_; private: + mutable TICLCandidateValidator candidateVal; CaloParticleSelector cpSelector; std::shared_ptr tools_; std::map cummatbudg; diff --git a/Validation/HGCalValidation/interface/TICLCandidateValidator.h b/Validation/HGCalValidation/interface/TICLCandidateValidator.h new file mode 100644 index 0000000000000..33713dd2e188a --- /dev/null +++ b/Validation/HGCalValidation/interface/TICLCandidateValidator.h @@ -0,0 +1,118 @@ +#ifndef Validation_HGCalValidation_TICLCandidateValidator_h +#define Validation_HGCalValidation_TICLCandidateValidator_h + +#include +#include +#include + +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" + +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterHitLCAssociator.h" + +#include "DQMServices/Core/interface/DQMStore.h" + +class TICLCandidateValidator { +public: + typedef dqm::legacy::DQMStore DQMStore; + typedef dqm::legacy::MonitorElement MonitorElement; + + TICLCandidateValidator(){}; + TICLCandidateValidator(edm::EDGetTokenT> TICLCandidates, + edm::EDGetTokenT> simTICLCandidatesToken, + edm::EDGetTokenT> recoTracksToken, + edm::EDGetTokenT> trackstersToken, + edm::EDGetTokenT associatorMapRtSToken, + edm::EDGetTokenT associatorMapStRToken); + ~TICLCandidateValidator(); + + void bookCandidatesHistos(DQMStore::IBooker& ibook, std::string baseDir); + + void fillCandidateHistos(const edm::Event& event, edm::Handle simTrackstersCP_h); + +private: + dqm::reco::MonitorElement* h_tracksters_in_candidate; + dqm::reco::MonitorElement* h_candidate_raw_energy; + dqm::reco::MonitorElement* h_candidate_regressed_energy; + dqm::reco::MonitorElement* h_candidate_pT; + dqm::reco::MonitorElement* h_candidate_charge; + dqm::reco::MonitorElement* h_candidate_pdgId; + dqm::reco::MonitorElement* h_candidate_partType; + + std::vector h_den_chg_energy_candidate; + std::vector h_num_chg_energy_candidate_track; + std::vector h_num_chg_energy_candidate_pdgId; + std::vector h_num_chg_energy_candidate_energy; + std::vector h_den_chg_pt_candidate; + std::vector h_num_chg_pt_candidate_track; + std::vector h_num_chg_pt_candidate_pdgId; + std::vector h_num_chg_pt_candidate_energy; + std::vector h_den_chg_eta_candidate; + std::vector h_num_chg_eta_candidate_track; + std::vector h_num_chg_eta_candidate_pdgId; + std::vector h_num_chg_eta_candidate_energy; + std::vector h_den_chg_phi_candidate; + std::vector h_num_chg_phi_candidate_track; + std::vector h_num_chg_phi_candidate_pdgId; + std::vector h_num_chg_phi_candidate_energy; + + std::vector h_den_neut_energy_candidate; + std::vector h_num_neut_energy_candidate_pdgId; + std::vector h_num_neut_energy_candidate_energy; + std::vector h_den_neut_pt_candidate; + std::vector h_num_neut_pt_candidate_pdgId; + std::vector h_num_neut_pt_candidate_energy; + std::vector h_den_neut_eta_candidate; + std::vector h_num_neut_eta_candidate_pdgId; + std::vector h_num_neut_eta_candidate_energy; + std::vector h_den_neut_phi_candidate; + std::vector h_num_neut_phi_candidate_pdgId; + std::vector h_num_neut_phi_candidate_energy; + + std::vector h_den_fake_chg_energy_candidate; + std::vector h_num_fake_chg_energy_candidate_track; + std::vector h_num_fake_chg_energy_candidate_pdgId; + std::vector h_num_fake_chg_energy_candidate_energy; + std::vector h_den_fake_chg_pt_candidate; + std::vector h_num_fake_chg_pt_candidate_track; + std::vector h_num_fake_chg_pt_candidate_pdgId; + std::vector h_num_fake_chg_pt_candidate_energy; + std::vector h_den_fake_chg_eta_candidate; + std::vector h_num_fake_chg_eta_candidate_track; + std::vector h_num_fake_chg_eta_candidate_pdgId; + std::vector h_num_fake_chg_eta_candidate_energy; + std::vector h_den_fake_chg_phi_candidate; + std::vector h_num_fake_chg_phi_candidate_track; + std::vector h_num_fake_chg_phi_candidate_pdgId; + std::vector h_num_fake_chg_phi_candidate_energy; + + std::vector h_den_fake_neut_energy_candidate; + std::vector h_num_fake_neut_energy_candidate_pdgId; + std::vector h_num_fake_neut_energy_candidate_energy; + std::vector h_den_fake_neut_pt_candidate; + std::vector h_num_fake_neut_pt_candidate_pdgId; + std::vector h_num_fake_neut_pt_candidate_energy; + std::vector h_den_fake_neut_eta_candidate; + std::vector h_num_fake_neut_eta_candidate_pdgId; + std::vector h_num_fake_neut_eta_candidate_energy; + std::vector h_den_fake_neut_phi_candidate; + std::vector h_num_fake_neut_phi_candidate_pdgId; + std::vector h_num_fake_neut_phi_candidate_energy; + + edm::EDGetTokenT> TICLCandidatesToken_; + edm::EDGetTokenT> simTICLCandidatesToken_; + edm::EDGetTokenT> recoTracksToken_; + edm::EDGetTokenT> trackstersToken_; + edm::EDGetTokenT associatorMapRtSToken_; + edm::EDGetTokenT associatorMapStRToken_; +}; + +#endif diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 3750312898480..9cd9dc27ebdc6 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -30,6 +30,8 @@ HGCalValidator::HGCalValidator(const edm::ParameterSet& pset) label_TSToCPLinking_(pset.getParameter("label_TSToCPLinking")), label_TSToSTSPR_(pset.getParameter("label_TSToSTSPR")), label_clustersmask(pset.getParameter>("LayerClustersInputMask")), + doCandidatesPlots_(pset.getUntrackedParameter("doCandidatesPlots")), + label_candidates_(pset.getParameter("ticlCandidates")), cummatbudinxo_(pset.getParameter("cummatbudinxo")) { //In this way we can easily generalize to associations between other objects also. const edm::InputTag& label_cp_effic_tag = pset.getParameter("label_cp_effic"); @@ -55,6 +57,29 @@ HGCalValidator::HGCalValidator(const edm::ParameterSet& pset) layerclusters_ = consumes(label_lcl); + if (doCandidatesPlots_) { + edm::EDGetTokenT> TICLCandidatesToken = + consumes>(pset.getParameter("ticlTrackstersMerge")); + edm::EDGetTokenT> simTICLCandidatesToken = + consumes>(pset.getParameter("simTiclCandidates")); + edm::EDGetTokenT> recoTracksToken = + consumes>(pset.getParameter("recoTracks")); + edm::EDGetTokenT> trackstersToken = + consumes>(pset.getParameter("ticlTrackstersMerge")); + //consumes>(pset.getParameter("trackstersclue3d")); + edm::EDGetTokenT associatorMapRtSToken = + consumes(pset.getParameter("mergeRecoToSimAssociator")); + edm::EDGetTokenT associatorMapStRToken = + consumes(pset.getParameter("mergeSimToRecoAssociator")); + + candidateVal = TICLCandidateValidator(TICLCandidatesToken, + simTICLCandidatesToken, + recoTracksToken, + trackstersToken, + associatorMapRtSToken, + associatorMapStRToken); + } + for (auto& itag : label_tst) { label_tstTokens.push_back(consumes(itag)); } @@ -219,6 +244,13 @@ void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook, ibook, histograms.histoProducerAlgo, HGVHistoProducerAlgo::validationType::PatternRecognition); } } //end of booking Tracksters loop + + // Booking histograms concerning TICL candidates + if (doCandidatesPlots_) { + ibook.cd(); + ibook.setCurrentFolder(dirName_ + label_candidates_.label()); + candidateVal.bookCandidatesHistos(ibook, dirName_ + label_candidates_.label()); + } } void HGCalValidator::cpParametersAndSelection(const Histograms& histograms, @@ -432,4 +464,9 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, totallayers_to_monitor_); } } //end of loop over Trackster input labels + + // tracksters histograms + if (doCandidatesPlots_) { + candidateVal.fillCandidateHistos(event, simTracksterFromCPHandle); + } } diff --git a/Validation/HGCalValidation/python/HGCalPostProcessor_cff.py b/Validation/HGCalValidation/python/HGCalPostProcessor_cff.py index d87f588a30a75..219f9ce6fefde 100644 --- a/Validation/HGCalValidation/python/HGCalPostProcessor_cff.py +++ b/Validation/HGCalValidation/python/HGCalPostProcessor_cff.py @@ -3,7 +3,7 @@ from Validation.HGCalValidation.HGCalSimHitsClient_cff import * from Validation.HGCalValidation.HGCalDigiClient_cff import * from Validation.HGCalValidation.HGCalRecHitsClient_cff import * -from Validation.HGCalValidation.PostProcessorHGCAL_cfi import postProcessorHGCALlayerclusters,postProcessorHGCALsimclusters,postProcessorHGCALTracksters +from Validation.HGCalValidation.PostProcessorHGCAL_cfi import postProcessorHGCALlayerclusters,postProcessorHGCALsimclusters,postProcessorHGCALTracksters,postProcessorHGCALCandidates hgcalPostProcessor = cms.Sequence(hgcalSimHitClientEE + hgcalSimHitClientHEF @@ -18,4 +18,5 @@ hgcalValidatorPostProcessor = cms.Sequence( postProcessorHGCALlayerclusters+ postProcessorHGCALsimclusters+ - postProcessorHGCALTracksters) + postProcessorHGCALTracksters+ + postProcessorHGCALCandidates) diff --git a/Validation/HGCalValidation/python/HGCalValidator_cfi.py b/Validation/HGCalValidation/python/HGCalValidator_cfi.py index db5ce31023bff..eae93bf63017e 100644 --- a/Validation/HGCalValidation/python/HGCalValidator_cfi.py +++ b/Validation/HGCalValidation/python/HGCalValidator_cfi.py @@ -51,6 +51,16 @@ label_TS = cms.string("Morphology"), label_TSToCPLinking = cms.string("TSToCP_linking"), label_TSToSTSPR = cms.string("TSToSTS_patternRecognition"), + #candidates plots + doCandidatesPlots = cms.untracked.bool(True), + ticlCandidates = cms.string("ticlCandidates"), + + ticlTrackstersMerge = cms.InputTag("ticlTrackstersMerge"), + simTiclCandidates = cms.InputTag("ticlSimTracksters"), + recoTracks = cms.InputTag("generalTracks"), + trackstersclue3d = cms.InputTag("ticlTrackstersCLUE3DHigh"), + mergeRecoToSimAssociator = cms.InputTag("tracksterSimTracksterAssociationLinking", "recoToSim"), + mergeSimToRecoAssociator = cms.InputTag("tracksterSimTracksterAssociationLinking", "simToReco"), #The cumulative material budget in front of each layer. To be more specific, it #is the material budget just in front of the active material (not including it). diff --git a/Validation/HGCalValidation/python/PostProcessorHGCAL_cfi.py b/Validation/HGCalValidation/python/PostProcessorHGCAL_cfi.py index 89b639d95739c..72f1b3de469d7 100644 --- a/Validation/HGCalValidation/python/PostProcessorHGCAL_cfi.py +++ b/Validation/HGCalValidation/python/PostProcessorHGCAL_cfi.py @@ -82,3 +82,39 @@ outputFileName = cms.untracked.string(""), verbose = cms.untracked.uint32(4) ) + +neutrals = ["photons", "neutral_pions", "neutral_hadrons"] +charged = ["electrons", "muons", "charged_hadrons"] +variables = ["energy", "pt", "eta", "phi"] +subDirsCandidates = [prefix + hgcalValidator.ticlCandidates.value() + "/" + c for cands in (neutrals, charged) for c in cands] +eff_candidates = [] + +for c in charged: + for var in variables: + eff_candidates.append("eff_"+c+"_track_"+var+" '"+c.replace("_", " ")+" candidates track efficiency vs "+var+"' num_track_cand_vs_"+var+"_"+c+" den_cand_vs_"+var+"_"+c) + eff_candidates.append("eff_"+c+"_pid_"+var+" '"+c.replace("_", " ")+" candidates track + pid efficiency vs "+var+"' num_pid_cand_vs_"+var+"_"+c+" den_cand_vs_"+var+"_"+c) + eff_candidates.append("eff_"+c+"_energy_"+var+" '"+c.replace("_", " ")+" candidates track + pid + energy efficiency vs "+var+"' num_energy_cand_vs_"+var+"_"+c+" den_cand_vs_"+var+"_"+c) +for n in neutrals: + for var in variables: + eff_candidates.append("eff_"+n+"_pid_"+var+" '"+n.replace("_", " ")+" candidates pid efficiency vs "+var+"' num_pid_cand_vs_"+var+"_"+n+" den_cand_vs_"+var+"_"+n) + eff_candidates.append("eff_"+n+"_energy_"+var+" '"+n.replace("_", " ")+" candidates pid + energy efficiency vs "+var+"' num_energy_cand_vs_"+var+"_"+n+" den_cand_vs_"+var+"_"+n) + +for c in charged: + for var in variables: + eff_candidates.append("fake_"+c+"_track_"+var+" '"+c.replace("_", " ")+" candidates track fake vs "+var+"' num_fake_track_cand_vs_"+var+"_"+c+" den_fake_cand_vs_"+var+"_"+c) + eff_candidates.append("fake_"+c+"_pid_"+var+" '"+c.replace("_", " ")+" candidates track + pid fake vs "+var+"' num_fake_pid_cand_vs_"+var+"_"+c+" den_fake_cand_vs_"+var+"_"+c) + eff_candidates.append("fake_"+c+"_energy_"+var+" '"+c.replace("_", " ")+" candidates track + pid + energy fake vs "+var+"' num_fake_energy_cand_vs_"+var+"_"+c+" den_fake_cand_vs_"+var+"_"+c) +for n in neutrals: + for var in variables: + eff_candidates.append("fake_"+n+"_pid_"+var+" '"+n.replace("_", " ")+" candidates pid fake vs "+var+"' num_fake_pid_cand_vs_"+var+"_"+n+" den_fake_cand_vs_"+var+"_"+n) + eff_candidates.append("fake_"+n+"_energy_"+var+" '"+n.replace("_", " ")+" candidates pid + energy fake vs "+var+"' num_fake_energy_cand_vs_"+var+"_"+n+" den_fake_cand_vs_"+var+"_"+n) + +postProcessorHGCALCandidates = DQMEDHarvester('DQMGenericClient', + subDirs = cms.untracked.vstring(subDirsCandidates), + efficiency = cms.vstring(eff_candidates), + resolution = cms.vstring(), + cumulativeDists = cms.untracked.vstring(), + noFlowDists = cms.untracked.vstring(), + outputFileName = cms.untracked.string(""), + verbose = cms.untracked.uint32(4) +) diff --git a/Validation/HGCalValidation/src/TICLCandidateValidator.cc b/Validation/HGCalValidation/src/TICLCandidateValidator.cc new file mode 100644 index 0000000000000..f69e9f4ada559 --- /dev/null +++ b/Validation/HGCalValidation/src/TICLCandidateValidator.cc @@ -0,0 +1,655 @@ +#include +#include +#include + +#include "Validation/HGCalValidation/interface/TICLCandidateValidator.h" +#include "RecoHGCal/TICL/interface/commons.h" + +TICLCandidateValidator::TICLCandidateValidator( + edm::EDGetTokenT> ticlCandidates, + edm::EDGetTokenT> simTICLCandidatesToken, + edm::EDGetTokenT> recoTracksToken, + edm::EDGetTokenT> trackstersToken, + edm::EDGetTokenT associatorMapRtSToken, + edm::EDGetTokenT associatorMapStRToken) + : TICLCandidatesToken_(ticlCandidates), + simTICLCandidatesToken_(simTICLCandidatesToken), + recoTracksToken_(recoTracksToken), + trackstersToken_(trackstersToken), + associatorMapRtSToken_(associatorMapRtSToken), + associatorMapStRToken_(associatorMapStRToken) {} + +TICLCandidateValidator::~TICLCandidateValidator() {} + +void TICLCandidateValidator::bookCandidatesHistos(DQMStore::IBooker& ibook, std::string baseDir) { + // book CAND histos + h_tracksters_in_candidate = ibook.book1D("N of tracksters in candidate", "N of tracksters in candidate", 100, 0, 99); + h_candidate_raw_energy = ibook.book1D("Candidates raw energy", "Candidates raw energy;E (GeV)", 250, 0, 250); + h_candidate_regressed_energy = + ibook.book1D("Candidates regressed energy", "Candidates regressed energy;E (GeV)", 250, 0, 250); + h_candidate_pT = ibook.book1D("Candidates pT", "Candidates pT;p_{T}", 250, 0, 250); + h_candidate_charge = ibook.book1D("Candidates charge", "Candidates charge;Charge", 3, -1.5, 1.5); + h_candidate_pdgId = ibook.book1D("Candidates PDG Id", "Candidates PDG ID", 100, -220, 220); + h_candidate_partType = ibook.book1D("Candidates type", "Candidates type", 9, -0.5, 8.5); + + // neutral: photon, pion, hadron + const std::vector neutrals{"photons", "neutral_pions", "neutral_hadrons"}; + for (long unsigned int i = 0; i < neutrals.size(); i++) { + ibook.setCurrentFolder(baseDir + "/" + neutrals[i]); + + h_den_fake_neut_energy_candidate.push_back( + ibook.book1D("den_fake_cand_vs_energy_" + neutrals[i], neutrals[i] + " vs energy;E (GeV)", 250, 0, 250)); + h_num_fake_neut_energy_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_energy_" + neutrals[i], neutrals[i] + " PID fake vs energy;E (GeV)", 250, 0, 250)); + h_num_fake_neut_energy_candidate_energy.push_back( + ibook.book1D("num_fake_energy_cand_vs_energy_" + neutrals[i], + neutrals[i] + " PID and energy fake vs energy;E (GeV)", + 250, + 0, + 250)); + h_den_fake_neut_pt_candidate.push_back( + ibook.book1D("den_fake_cand_vs_pt_" + neutrals[i], neutrals[i] + " vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_fake_neut_pt_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " PID fake vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_fake_neut_pt_candidate_energy.push_back(ibook.book1D("num_fake_energy_cand_vs_pt_" + neutrals[i], + neutrals[i] + " PID and energy fake vs pT;p_{T} (GeV)", + 250, + 0, + 250)); + h_den_fake_neut_eta_candidate.push_back( + ibook.book1D("den_fake_cand_vs_eta_" + neutrals[i], neutrals[i] + " vs eta;#eta (GeV)", 100, -3, 3)); + h_num_fake_neut_eta_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " PID fake vs eta;#eta (GeV)", 100, -3, 3)); + h_num_fake_neut_eta_candidate_energy.push_back(ibook.book1D("num_fake_energy_cand_vs_eta_" + neutrals[i], + neutrals[i] + " PID and energy fake vs eta;#eta (GeV)", + 100, + -3, + 3)); + h_den_fake_neut_phi_candidate.push_back(ibook.book1D( + "den_fake_cand_vs_phi_" + neutrals[i], neutrals[i] + " vs phi;#phi (GeV)", 100, -3.14159, 3.14159)); + h_num_fake_neut_phi_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_phi_" + neutrals[i], neutrals[i] + " PID fake vs phi;#phi (GeV)", 100, -3.14159, 3.14159)); + h_num_fake_neut_phi_candidate_energy.push_back(ibook.book1D("num_fake_energy_cand_vs_phi_" + neutrals[i], + neutrals[i] + " PID and energy fake vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + + h_den_neut_energy_candidate.push_back( + ibook.book1D("den_cand_vs_energy_" + neutrals[i], neutrals[i] + " vs energy;E (GeV)", 250, 0, 250)); + h_num_neut_energy_candidate_pdgId.push_back( + ibook.book1D("num_pid_cand_vs_energy_" + neutrals[i], + neutrals[i] + " track and PID efficiency vs energy;E (GeV)", + 250, + 0, + 250)); + h_num_neut_energy_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_energy_" + neutrals[i], + neutrals[i] + " track, PID and energy efficiency vs energy;E (GeV)", + 250, + 0, + 250)); + h_den_neut_pt_candidate.push_back( + ibook.book1D("den_cand_vs_pt_" + neutrals[i], neutrals[i] + " vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_neut_pt_candidate_pdgId.push_back(ibook.book1D( + "num_pid_cand_vs_pt_" + neutrals[i], neutrals[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_neut_pt_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_pt_" + neutrals[i], + neutrals[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)", + 250, + 0, + 250)); + h_den_neut_eta_candidate.push_back( + ibook.book1D("den_cand_vs_eta_" + neutrals[i], neutrals[i] + " vs eta;#eta (GeV)", 100, -3, 3)); + h_num_neut_eta_candidate_pdgId.push_back(ibook.book1D( + "num_pid_cand_vs_eta_" + neutrals[i], neutrals[i] + " track and PID efficiency vs eta;#eta (GeV)", 100, -3, 3)); + h_num_neut_eta_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_eta_" + neutrals[i], + neutrals[i] + " track, PID and energy efficiency vs eta;#eta (GeV)", + 100, + -3, + 3)); + h_den_neut_phi_candidate.push_back( + ibook.book1D("den_cand_vs_phi_" + neutrals[i], neutrals[i] + " vs phi;#phi (GeV)", 100, -3.14159, 3.14159)); + h_num_neut_phi_candidate_pdgId.push_back(ibook.book1D("num_pid_cand_vs_phi_" + neutrals[i], + neutrals[i] + " track and PID efficiency vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + h_num_neut_phi_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_phi_" + neutrals[i], + neutrals[i] + " track, PID and energy efficiency vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + } + // charged: electron, muon, hadron + const std::vector charged{"electrons", "muons", "charged_hadrons"}; + for (long unsigned int i = 0; i < charged.size(); i++) { + ibook.setCurrentFolder(baseDir + "/" + charged[i]); + + h_den_fake_chg_energy_candidate.push_back( + ibook.book1D("den_fake_cand_vs_energy_" + charged[i], charged[i] + " vs energy;E (GeV)", 250, 0, 250)); + h_num_fake_chg_energy_candidate_track.push_back(ibook.book1D( + "num_fake_track_cand_vs_energy_" + charged[i], charged[i] + " track fake vs energy;E (GeV)", 250, 0, 250)); + h_num_fake_chg_energy_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_energy_" + charged[i], charged[i] + " track and PID fake vs energy;E (GeV)", 250, 0, 250)); + h_num_fake_chg_energy_candidate_energy.push_back( + ibook.book1D("num_fake_energy_cand_vs_energy_" + charged[i], + charged[i] + " track, PID and energy fake vs energy;E (GeV)", + 250, + 0, + 250)); + h_den_fake_chg_pt_candidate.push_back( + ibook.book1D("den_fake_cand_vs_pt_" + charged[i], charged[i] + " vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_fake_chg_pt_candidate_track.push_back(ibook.book1D( + "num_fake_track_cand_vs_pt_" + charged[i], charged[i] + " track fake vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_fake_chg_pt_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID fake vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_fake_chg_pt_candidate_energy.push_back( + ibook.book1D("num_fake_energy_cand_vs_pt_" + charged[i], + charged[i] + " track, PID and energy fake vs pT;p_{T} (GeV)", + 250, + 0, + 250)); + h_den_fake_chg_eta_candidate.push_back( + ibook.book1D("den_fake_cand_vs_eta_" + charged[i], charged[i] + " vs eta;#eta (GeV)", 100, -3, 3)); + h_num_fake_chg_eta_candidate_track.push_back(ibook.book1D( + "num_fake_track_cand_vs_eta_" + charged[i], charged[i] + " track fake vs eta;#eta (GeV)", 100, -3, 3)); + h_num_fake_chg_eta_candidate_pdgId.push_back(ibook.book1D( + "num_fake_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID fake vs eta;#eta (GeV)", 100, -3, 3)); + h_num_fake_chg_eta_candidate_energy.push_back( + ibook.book1D("num_fake_energy_cand_vs_eta_" + charged[i], + charged[i] + " track, PID and energy fake vs eta;#eta (GeV)", + 100, + -3, + 3)); + h_den_fake_chg_phi_candidate.push_back( + ibook.book1D("den_fake_cand_vs_phi_" + charged[i], charged[i] + " vs phi;#phi (GeV)", 100, -3.14159, 3.14159)); + h_num_fake_chg_phi_candidate_track.push_back(ibook.book1D("num_fake_track_cand_vs_phi_" + charged[i], + charged[i] + " track fake vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + h_num_fake_chg_phi_candidate_pdgId.push_back(ibook.book1D("num_fake_pid_cand_vs_phi_" + charged[i], + charged[i] + " track and PID fake vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + h_num_fake_chg_phi_candidate_energy.push_back( + ibook.book1D("num_fake_energy_cand_vs_phi_" + charged[i], + charged[i] + " track, PID and energy fake vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + + h_den_chg_energy_candidate.push_back( + ibook.book1D("den_cand_vs_energy_" + charged[i], charged[i] + " vs energy;E (GeV)", 250, 0, 250)); + h_num_chg_energy_candidate_track.push_back(ibook.book1D( + "num_track_cand_vs_energy_" + charged[i], charged[i] + " track efficiency vs energy;E (GeV)", 250, 0, 250)); + h_num_chg_energy_candidate_pdgId.push_back(ibook.book1D("num_pid_cand_vs_energy_" + charged[i], + charged[i] + " track and PID efficiency vs energy;E (GeV)", + 250, + 0, + 250)); + h_num_chg_energy_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_energy_" + charged[i], + charged[i] + " track, PID and energy efficiency vs energy;E (GeV)", + 250, + 0, + 250)); + h_den_chg_pt_candidate.push_back( + ibook.book1D("den_cand_vs_pt_" + charged[i], charged[i] + " vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_chg_pt_candidate_track.push_back(ibook.book1D( + "num_track_cand_vs_pt_" + charged[i], charged[i] + " track efficiency vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_chg_pt_candidate_pdgId.push_back(ibook.book1D( + "num_pid_cand_vs_pt_" + charged[i], charged[i] + " track and PID efficiency vs pT;p_{T} (GeV)", 250, 0, 250)); + h_num_chg_pt_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_pt_" + charged[i], + charged[i] + " track, PID and energy efficiency vs pT;p_{T} (GeV)", + 250, + 0, + 250)); + h_den_chg_eta_candidate.push_back( + ibook.book1D("den_cand_vs_eta_" + charged[i], charged[i] + " vs eta;#eta (GeV)", 100, -3, 3)); + h_num_chg_eta_candidate_track.push_back(ibook.book1D( + "num_track_cand_vs_eta_" + charged[i], charged[i] + " track efficiency vs eta;#eta (GeV)", 100, -3, 3)); + h_num_chg_eta_candidate_pdgId.push_back(ibook.book1D( + "num_pid_cand_vs_eta_" + charged[i], charged[i] + " track and PID efficiency vs eta;#eta (GeV)", 100, -3, 3)); + h_num_chg_eta_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_eta_" + charged[i], + charged[i] + " track, PID and energy efficiency vs eta;#eta (GeV)", + 100, + -3, + 3)); + h_den_chg_phi_candidate.push_back( + ibook.book1D("den_cand_vs_phi_" + charged[i], charged[i] + " vs phi;#phi (GeV)", 100, -3.14159, 3.14159)); + h_num_chg_phi_candidate_track.push_back(ibook.book1D("num_track_cand_vs_phi_" + charged[i], + charged[i] + " track efficiency vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + h_num_chg_phi_candidate_pdgId.push_back(ibook.book1D("num_pid_cand_vs_phi_" + charged[i], + charged[i] + " track and PID efficiency vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + h_num_chg_phi_candidate_energy.push_back( + ibook.book1D("num_energy_cand_vs_phi_" + charged[i], + charged[i] + " track, PID and energy efficiency vs phi;#phi (GeV)", + 100, + -3.14159, + 3.14159)); + } +} + +void TICLCandidateValidator::fillCandidateHistos(const edm::Event& event, + edm::Handle simTrackstersCP_h) { + edm::Handle> TICLCandidates_h; + event.getByToken(TICLCandidatesToken_, TICLCandidates_h); + auto TICLCandidates = *TICLCandidates_h; + + edm::Handle> simTICLCandidates_h; + event.getByToken(simTICLCandidatesToken_, simTICLCandidates_h); + auto simTICLCandidates = *simTICLCandidates_h; + + edm::Handle> recoTracks_h; + event.getByToken(recoTracksToken_, recoTracks_h); + auto recoTracks = *recoTracks_h; + + edm::Handle> Tracksters_h; + event.getByToken(trackstersToken_, Tracksters_h); + auto trackstersMerged = *Tracksters_h; + + edm::Handle mergeTsRecoToSim_h; + event.getByToken(associatorMapRtSToken_, mergeTsRecoToSim_h); + auto const& mergeTsRecoToSimMap = *mergeTsRecoToSim_h; + + edm::Handle mergeTsSimToReco_h; + event.getByToken(associatorMapStRToken_, mergeTsSimToReco_h); + auto const& mergeTsSimToRecoMap = *mergeTsSimToReco_h; + + // candidates plots + for (const auto& cand : TICLCandidates) { + h_tracksters_in_candidate->Fill(cand.tracksters().size()); + h_candidate_raw_energy->Fill(cand.rawEnergy()); + h_candidate_regressed_energy->Fill(cand.energy()); + h_candidate_pT->Fill(cand.pt()); + h_candidate_charge->Fill(cand.charge()); + h_candidate_pdgId->Fill(cand.pdgId()); + const auto& arr = cand.idProbabilities(); + h_candidate_partType->Fill(std::max_element(arr.begin(), arr.end()) - arr.begin()); + } + + std::vector chargedCandidates; + std::vector neutralCandidates; + chargedCandidates.reserve(simTICLCandidates.size()); + neutralCandidates.reserve(simTICLCandidates.size()); + + for (size_t i = 0; i < simTICLCandidates.size(); ++i) { + const auto& simCand = simTICLCandidates[i]; + const auto particleType = ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), simCand.charge()); + if (particleType == ticl::Trackster::ParticleType::electron or + particleType == ticl::Trackster::ParticleType::muon or + particleType == ticl::Trackster::ParticleType::charged_hadron) + chargedCandidates.emplace_back(i); + else if (particleType == ticl::Trackster::ParticleType::photon or + particleType == ticl::Trackster::ParticleType::neutral_pion or + particleType == ticl::Trackster::ParticleType::neutral_hadron) + neutralCandidates.emplace_back(i); + // should consider also unknown ? + } + + chargedCandidates.shrink_to_fit(); + neutralCandidates.shrink_to_fit(); + + for (const auto i : chargedCandidates) { + const auto& simCand = simTICLCandidates[i]; + auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 1))); + /* 11 (type 1) becomes 0 + * 13 (type 2) becomes 1 + * 211 (type 4) becomes 2 + */ + int32_t simCandTrackIdx = -1; + if (simCand.trackPtr().get() != nullptr) + simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr(recoTracks_h, 0).get(); + else { + // no reco track, but simCand is charged + continue; + } + if (simCand.trackPtr().get()->pt() < 1 or simCand.trackPtr().get()->missingOuterHits() > 5 or + not simCand.trackPtr().get()->quality(reco::TrackBase::highPurity)) + continue; + + // +1 to all denominators + h_den_chg_energy_candidate[index]->Fill(simCand.rawEnergy()); + h_den_chg_pt_candidate[index]->Fill(simCand.pt()); + h_den_chg_eta_candidate[index]->Fill(simCand.eta()); + h_den_chg_phi_candidate[index]->Fill(simCand.phi()); + + int32_t cand_idx = -1; + const edm::Ref stsRef(simTrackstersCP_h, i); + const auto ts_iter = mergeTsSimToRecoMap.find(stsRef); + float shared_energy = 0.; + // search for reco cand associated + if (ts_iter != mergeTsSimToRecoMap.end()) { + const auto& tsAssoc = (ts_iter->val); + std::vector MergeTracksters_simToReco; + std::vector MergeTracksters_simToReco_score; + std::vector MergeTracksters_simToReco_sharedE; + MergeTracksters_simToReco.reserve(tsAssoc.size()); + MergeTracksters_simToReco_score.reserve(tsAssoc.size()); + MergeTracksters_simToReco_sharedE.reserve(tsAssoc.size()); + for (auto& ts : tsAssoc) { + auto ts_id = (ts.first).get() - (edm::Ref(Tracksters_h, 0)).get(); + MergeTracksters_simToReco.push_back(ts_id); + MergeTracksters_simToReco_score.push_back(ts.second.second); + MergeTracksters_simToReco_sharedE.push_back(ts.second.first); + } + auto min_idx = std::min_element(MergeTracksters_simToReco_score.begin(), MergeTracksters_simToReco_score.end()); + if (*min_idx != 1) { + cand_idx = MergeTracksters_simToReco[min_idx - MergeTracksters_simToReco_score.begin()]; + shared_energy = MergeTracksters_simToReco_sharedE[min_idx - MergeTracksters_simToReco_score.begin()]; + } + } + + // no reco associated to sim + if (cand_idx == -1) + continue; + + const auto& recoCand = TICLCandidates[cand_idx]; + if (recoCand.trackPtr().get() != nullptr) { + const auto candTrackIdx = recoCand.trackPtr().get() - edm::Ptr(recoTracks_h, 0).get(); + if (simCandTrackIdx == candTrackIdx) { + // +1 to track num + h_num_chg_energy_candidate_track[index]->Fill(simCand.rawEnergy()); + h_num_chg_pt_candidate_track[index]->Fill(simCand.pt()); + h_num_chg_eta_candidate_track[index]->Fill(simCand.eta()); + h_num_chg_phi_candidate_track[index]->Fill(simCand.phi()); + } else { + continue; + } + } else { + continue; + } + + //step 2: PID + if (simCand.pdgId() == recoCand.pdgId()) { + // +1 to num pdg id + h_num_chg_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy()); + h_num_chg_pt_candidate_pdgId[index]->Fill(simCand.pt()); + h_num_chg_eta_candidate_pdgId[index]->Fill(simCand.eta()); + h_num_chg_phi_candidate_pdgId[index]->Fill(simCand.phi()); + + //step 3: energy + if (shared_energy / simCand.rawEnergy() > 0.5) { + // +1 to ene num + h_num_chg_energy_candidate_energy[index]->Fill(simCand.rawEnergy()); + h_num_chg_pt_candidate_energy[index]->Fill(simCand.pt()); + h_num_chg_eta_candidate_energy[index]->Fill(simCand.eta()); + h_num_chg_phi_candidate_energy[index]->Fill(simCand.phi()); + } + } + } + + for (const auto i : neutralCandidates) { + const auto& simCand = simTICLCandidates[i]; + auto index = int(ticl::tracksterParticleTypeFromPdgId(simCand.pdgId(), 0)) / 2; + /* 22 (type 0) becomes 0 + * 111 (type 3) becomes 1 + * 130 (type 5) becomes 2 + */ + //if (simCand.trackPtr().get() != nullptr) + // std::cout << "ERROR: NEUTRAL WITH TRACK\n"; + + h_den_neut_energy_candidate[index]->Fill(simCand.rawEnergy()); + h_den_neut_pt_candidate[index]->Fill(simCand.pt()); + h_den_neut_eta_candidate[index]->Fill(simCand.eta()); + h_den_neut_phi_candidate[index]->Fill(simCand.phi()); + + int32_t cand_idx = -1; + const edm::Ref stsRef(simTrackstersCP_h, i); + const auto ts_iter = mergeTsSimToRecoMap.find(stsRef); + float shared_energy = 0.; + // search for reco cand associated + if (ts_iter != mergeTsSimToRecoMap.end()) { + const auto& tsAssoc = (ts_iter->val); + std::vector MergeTracksters_simToReco; + std::vector MergeTracksters_simToReco_score; + std::vector MergeTracksters_simToReco_sharedE; + MergeTracksters_simToReco.reserve(tsAssoc.size()); + MergeTracksters_simToReco_score.reserve(tsAssoc.size()); + MergeTracksters_simToReco_sharedE.reserve(tsAssoc.size()); + for (auto& ts : tsAssoc) { + auto ts_id = (ts.first).get() - (edm::Ref(Tracksters_h, 0)).get(); + MergeTracksters_simToReco.push_back(ts_id); + MergeTracksters_simToReco_score.push_back(ts.second.second); + MergeTracksters_simToReco_sharedE.push_back(ts.second.first); + } + auto min_idx = std::min_element(MergeTracksters_simToReco_score.begin(), MergeTracksters_simToReco_score.end()); + if (*min_idx != 1) { + cand_idx = MergeTracksters_simToReco[min_idx - MergeTracksters_simToReco_score.begin()]; + shared_energy = MergeTracksters_simToReco_sharedE[min_idx - MergeTracksters_simToReco_score.begin()]; + } + } + + // no reco associated to sim + if (cand_idx == -1) + continue; + + const auto& recoCand = TICLCandidates[cand_idx]; + if (recoCand.trackPtr().get() != nullptr) + continue; + + //step 2: PID + if (simCand.pdgId() == recoCand.pdgId()) { + // +1 to num pdg id + h_num_neut_energy_candidate_pdgId[index]->Fill(simCand.rawEnergy()); + h_num_neut_pt_candidate_pdgId[index]->Fill(simCand.pt()); + h_num_neut_eta_candidate_pdgId[index]->Fill(simCand.eta()); + h_num_neut_phi_candidate_pdgId[index]->Fill(simCand.phi()); + + //step 3: energy + if (shared_energy / simCand.rawEnergy() > 0.5) { + // +1 to ene num + h_num_neut_energy_candidate_energy[index]->Fill(simCand.rawEnergy()); + h_num_neut_pt_candidate_energy[index]->Fill(simCand.pt()); + h_num_neut_eta_candidate_energy[index]->Fill(simCand.eta()); + h_num_neut_phi_candidate_energy[index]->Fill(simCand.phi()); + } + } + } + + // FAKE rate + chargedCandidates.clear(); + neutralCandidates.clear(); + chargedCandidates.reserve(TICLCandidates.size()); + neutralCandidates.reserve(TICLCandidates.size()); + + auto isCharged = [](int pdgId) { + pdgId = std::abs(pdgId); + return (pdgId == 11 or pdgId == 211 or pdgId == 13); + }; + + for (size_t i = 0; i < TICLCandidates.size(); ++i) { + const auto& cand = TICLCandidates[i]; + // const auto& arr = cand.idProbabilities(); + // const auto particleType = static_cast(std::max_element(arr.begin(), arr.end()) - arr.begin()); + // if (particleType == ticl::Trackster::ParticleType::electron or particleType == ticl::Trackster::ParticleType::muon or + // particleType == ticl::Trackster::ParticleType::charged_hadron) + // chargedCandidates.emplace_back(i); + // else if (particleType == ticl::Trackster::ParticleType::photon or particleType == ticl::Trackster::ParticleType::neutral_pion or + // particleType == ticl::Trackster::ParticleType::neutral_hadron) + // neutralCandidates.emplace_back(i); + const auto& charged = isCharged(cand.pdgId()); + if (charged) + chargedCandidates.emplace_back(i); + else + neutralCandidates.emplace_back(i); + + // should consider also unknown ? + } + + chargedCandidates.shrink_to_fit(); + neutralCandidates.shrink_to_fit(); + + // loop on charged + for (const auto i : chargedCandidates) { + const auto& cand = TICLCandidates[i]; + // const auto& arr = cand.idProbabilities(); + // auto index = std::log2(int(std::max_element(arr.begin(), arr.end()) - arr.begin())); + auto index = std::log2(int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 1))); + /* 11 (type 1) becomes 0 + * 13 (type 2) becomes 1 + * 211 (type 4) becomes 2 + */ + int32_t candTrackIdx = -1; + if (cand.trackPtr().get() != nullptr) + candTrackIdx = cand.trackPtr().get() - edm::Ptr(recoTracks_h, 0).get(); + else + continue; + // +1 to all denominators + h_den_fake_chg_energy_candidate[index]->Fill(cand.rawEnergy()); + h_den_fake_chg_pt_candidate[index]->Fill(cand.pt()); + h_den_fake_chg_eta_candidate[index]->Fill(cand.eta()); + h_den_fake_chg_phi_candidate[index]->Fill(cand.phi()); + + int32_t simCand_idx = -1; + const edm::Ref tsRef(Tracksters_h, i); + const auto sts_iter = mergeTsRecoToSimMap.find(tsRef); + float shared_energy = 0.; + // search for reco cand associated + if (sts_iter != mergeTsRecoToSimMap.end()) { + const auto& stsAssoc = (sts_iter->val); + std::vector MergeTracksters_recoToSim; + std::vector MergeTracksters_recoToSim_score; + std::vector MergeTracksters_recoToSim_sharedE; + MergeTracksters_recoToSim.reserve(stsAssoc.size()); + MergeTracksters_recoToSim_score.reserve(stsAssoc.size()); + MergeTracksters_recoToSim_sharedE.reserve(stsAssoc.size()); + for (auto& sts : stsAssoc) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); + MergeTracksters_recoToSim.push_back(sts_id); + MergeTracksters_recoToSim_score.push_back(sts.second.second); + MergeTracksters_recoToSim_sharedE.push_back(sts.second.first); + } + auto min_idx = std::min_element(MergeTracksters_recoToSim_score.begin(), MergeTracksters_recoToSim_score.end()); + if (*min_idx != 1) { + simCand_idx = MergeTracksters_recoToSim[min_idx - MergeTracksters_recoToSim_score.begin()]; + shared_energy = MergeTracksters_recoToSim_sharedE[min_idx - MergeTracksters_recoToSim_score.begin()]; + } + } + + if (simCand_idx == -1) + continue; + + const auto& simCand = simTICLCandidates[simCand_idx]; + if (simCand.trackPtr().get() != nullptr) { + const auto simCandTrackIdx = simCand.trackPtr().get() - edm::Ptr(recoTracks_h, 0).get(); + if (simCandTrackIdx != candTrackIdx) { + // fake += 1 + h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy()); + h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt()); + h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta()); + h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi()); + continue; + } + } else { + // fake += 1 + h_num_fake_chg_energy_candidate_track[index]->Fill(cand.rawEnergy()); + h_num_fake_chg_pt_candidate_track[index]->Fill(cand.pt()); + h_num_fake_chg_eta_candidate_track[index]->Fill(cand.eta()); + h_num_fake_chg_phi_candidate_track[index]->Fill(cand.phi()); + continue; + } + + //step 2: PID + if (simCand.pdgId() != cand.pdgId()) { + // +1 to num fake pdg id + h_num_fake_chg_energy_candidate_pdgId[index]->Fill(cand.rawEnergy()); + h_num_fake_chg_pt_candidate_pdgId[index]->Fill(cand.pt()); + h_num_fake_chg_eta_candidate_pdgId[index]->Fill(cand.eta()); + h_num_fake_chg_phi_candidate_pdgId[index]->Fill(cand.phi()); + continue; + } + + //step 3: energy + if (shared_energy / simCand.rawEnergy() < 0.5) { + // +1 to ene num + h_num_fake_chg_energy_candidate_energy[index]->Fill(cand.rawEnergy()); + h_num_fake_chg_pt_candidate_energy[index]->Fill(cand.pt()); + h_num_fake_chg_eta_candidate_energy[index]->Fill(cand.eta()); + h_num_fake_chg_phi_candidate_energy[index]->Fill(cand.phi()); + } + } + // loop on neutrals + for (const auto i : neutralCandidates) { + const auto& cand = TICLCandidates[i]; + // const auto& arr = cand.idProbabilities(); + // auto index = int(std::max_element(arr.begin(), arr.end()) - arr.begin()) / 2; + auto index = int(ticl::tracksterParticleTypeFromPdgId(cand.pdgId(), 0)) / 2; + /* 22 (type 0) becomes 0 + * 111 (type 3) becomes 1 + * 130 (type 5) becomes 2 + */ + if (cand.trackPtr().get() != nullptr) + continue; + + // +1 to all denominators + h_den_fake_neut_energy_candidate[index]->Fill(cand.rawEnergy()); + h_den_fake_neut_pt_candidate[index]->Fill(cand.pt()); + h_den_fake_neut_eta_candidate[index]->Fill(cand.eta()); + h_den_fake_neut_phi_candidate[index]->Fill(cand.phi()); + + int32_t simCand_idx = -1; + const edm::Ref tsRef(Tracksters_h, i); + const auto sts_iter = mergeTsRecoToSimMap.find(tsRef); + float shared_energy = 0.; + // search for reco cand associated + if (sts_iter != mergeTsRecoToSimMap.end()) { + const auto& stsAssoc = (sts_iter->val); + std::vector MergeTracksters_recoToSim; + std::vector MergeTracksters_recoToSim_score; + std::vector MergeTracksters_recoToSim_sharedE; + MergeTracksters_recoToSim.reserve(stsAssoc.size()); + MergeTracksters_recoToSim_score.reserve(stsAssoc.size()); + MergeTracksters_recoToSim_sharedE.reserve(stsAssoc.size()); + for (auto& sts : stsAssoc) { + auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); + MergeTracksters_recoToSim.push_back(sts_id); + MergeTracksters_recoToSim_score.push_back(sts.second.second); + MergeTracksters_recoToSim_sharedE.push_back(sts.second.first); + } + auto min_idx = std::min_element(MergeTracksters_recoToSim_score.begin(), MergeTracksters_recoToSim_score.end()); + if (*min_idx != 1) { + simCand_idx = MergeTracksters_recoToSim[min_idx - MergeTracksters_recoToSim_score.begin()]; + shared_energy = MergeTracksters_recoToSim_sharedE[min_idx - MergeTracksters_recoToSim_score.begin()]; + } + } + + if (simCand_idx == -1) + continue; + + const auto& simCand = simTICLCandidates[simCand_idx]; + // check su simTICLCandidate track? + + //step 2: PID + if (simCand.pdgId() != cand.pdgId()) { + // +1 to num fake pdg id + h_num_fake_neut_energy_candidate_pdgId[index]->Fill(cand.rawEnergy()); + h_num_fake_neut_pt_candidate_pdgId[index]->Fill(cand.pt()); + h_num_fake_neut_eta_candidate_pdgId[index]->Fill(cand.eta()); + h_num_fake_neut_phi_candidate_pdgId[index]->Fill(cand.phi()); + continue; + } + + //step 3: energy + if (shared_energy / simCand.rawEnergy() < 0.5) { + // +1 to ene num + h_num_fake_neut_energy_candidate_energy[index]->Fill(cand.rawEnergy()); + h_num_fake_neut_pt_candidate_energy[index]->Fill(cand.pt()); + h_num_fake_neut_eta_candidate_energy[index]->Fill(cand.eta()); + h_num_fake_neut_phi_candidate_energy[index]->Fill(cand.phi()); + } + } +}