diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc index 6428196e2c096..ea0a56a852147 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc @@ -48,6 +48,9 @@ void fillOutputBranches(LSTEvent* event) { unsigned int n_accepted_simtrk = setSimTrackContainerBranches(event); + if (ana.jet_branches) + setGenJetBranches(event); + if (ana.occ_branches) setOccupancyBranches(event); @@ -192,12 +195,13 @@ void createT4DNNBranches() { //________________________________________________________________________________________________________________________________ void createJetBranches() { - ana.tx->createBranch>("sim_deltaEta"); - ana.tx->createBranch>("sim_deltaPhi"); - ana.tx->createBranch>("sim_deltaR"); - ana.tx->createBranch>("sim_jet_eta"); - ana.tx->createBranch>("sim_jet_phi"); - ana.tx->createBranch>("sim_jet_pt"); + ana.tx->createBranch>("sim_genjet_deltaEta"); + ana.tx->createBranch>("sim_genjet_deltaPhi"); + ana.tx->createBranch>("sim_genjet_deltaR"); + ana.tx->createBranch>("sim_genjet_idx"); + ana.tx->createBranch>("genjet_eta"); + ana.tx->createBranch>("genjet_phi"); + ana.tx->createBranch>("genjet_pt"); } //________________________________________________________________________________________________________________________________ @@ -616,6 +620,27 @@ void createOccupancyBranches() { ana.tx->createBranch("pT5_occupancies"); } +//________________________________________________________________________________________________________________________________ +void setGenJetBranches(LSTEvent* event) { + //-------------------------------------------- + // + // + // Gen Jets + // + // + //-------------------------------------------- + + auto const& trk_genjet_eta = trk.getVF("genjet_eta"); + auto const& trk_genjet_phi = trk.getVF("genjet_phi"); + auto const& trk_genjet_pt = trk.getVF("genjet_pt"); + + for (unsigned int ijet = 0; ijet < trk_genjet_pt.size(); ++ijet) { + ana.tx->pushbackToBranch("genjet_eta", trk_genjet_eta[ijet]); + ana.tx->pushbackToBranch("genjet_phi", trk_genjet_phi[ijet]); + ana.tx->pushbackToBranch("genjet_pt", trk_genjet_pt[ijet]); + } +} + //________________________________________________________________________________________________________________________________ unsigned int setSimTrackContainerBranches(LSTEvent* event) { //-------------------------------------------- @@ -669,19 +694,15 @@ unsigned int setSimTrackContainerBranches(LSTEvent* event) { // Now we have a list of "accepted" tracks (no condition on vtx_z/perp, nor pt, eta etc are applied yet) if (ana.jet_branches) { - auto const& trk_sim_deltaEta = trk.getVF("sim_deltaEta"); - auto const& trk_sim_deltaPhi = trk.getVF("sim_deltaPhi"); - auto const& trk_sim_deltaR = trk.getVF("sim_deltaR"); - auto const& trk_sim_jet_eta = trk.getVF("sim_jet_eta"); - auto const& trk_sim_jet_phi = trk.getVF("sim_jet_phi"); - auto const& trk_sim_jet_pt = trk.getVF("sim_jet_pt"); - - ana.tx->pushbackToBranch("sim_deltaEta", trk_sim_deltaEta[isimtrk]); - ana.tx->pushbackToBranch("sim_deltaPhi", trk_sim_deltaPhi[isimtrk]); - ana.tx->pushbackToBranch("sim_deltaR", trk_sim_deltaR[isimtrk]); - ana.tx->pushbackToBranch("sim_jet_eta", trk_sim_jet_eta[isimtrk]); - ana.tx->pushbackToBranch("sim_jet_phi", trk_sim_jet_phi[isimtrk]); - ana.tx->pushbackToBranch("sim_jet_pt", trk_sim_jet_pt[isimtrk]); + auto const& trk_sim_genjet_deltaEta = trk.getVF("sim_genjet_deltaEta"); + auto const& trk_sim_genjet_deltaPhi = trk.getVF("sim_genjet_deltaPhi"); + auto const& trk_sim_genjet_deltaR = trk.getVF("sim_genjet_deltaR"); + auto const& trk_sim_genjet_idx = trk.getVI("sim_genjet_idx"); + + ana.tx->pushbackToBranch("sim_genjet_deltaEta", trk_sim_genjet_deltaEta[isimtrk]); + ana.tx->pushbackToBranch("sim_genjet_deltaPhi", trk_sim_genjet_deltaPhi[isimtrk]); + ana.tx->pushbackToBranch("sim_genjet_deltaR", trk_sim_genjet_deltaR[isimtrk]); + ana.tx->pushbackToBranch("sim_genjet_idx", trk_sim_genjet_idx[isimtrk]); } // Fill the branch with simulated tracks. diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h index ee3786782200c..8a809576c080f 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h +++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h @@ -33,6 +33,7 @@ void createOccupancyBranches(); void fillOutputBranches(LSTEvent* event); void setOccupancyBranches(LSTEvent* event); +void setGenJetBranches(LSTEvent* event); unsigned int setSimTrackContainerBranches(LSTEvent* event); void setTrackCandidateBranches(LSTEvent* event, unsigned int n_accepted_tracks, diff --git a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py index d7b6aee0da42d..1d35f00af05f0 100755 --- a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py +++ b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py @@ -9,7 +9,7 @@ sel_choices = ["base", "loweta", "xtr", "vtr", "none"] metric_choices = ["eff", "fakerate", "duplrate", "fakeorduplrate", "avgOTlen"] -variable_choices = ["pt", "ptmtv", "ptlow", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "jet_eta", "jet_phi", "jet_pt"] +variable_choices = ["pt", "ptmtv", "ptlow", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "genJetEta", "genJetPhi", "genJetPt"] objecttype_choices = ["TC", "pT5", "T5", "pT3", "pLS", "T4", "pT5_lower", "pT3_lower", "T5_lower", "pLS_lower"] #lowerObjectType = ["pT5_lower", "pT3_lower", "T5_lower"] @@ -478,11 +478,11 @@ def set_label(eff, output_name, raw_number): title = "#phi diffs" elif "_deltaR" in output_name: title = "#Delta R" - elif "_jet_eta" in output_name: + elif "_genJetEta" in output_name: title = "jet #eta" - elif "_jet_phi" in output_name: + elif "_genJetPhi" in output_name: title = "jet #phi" - elif "_jet_pt" in output_name: + elif "_genJetPt" in output_name: title = "jet pT" elif "_dz" in output_name: title = "z [cm]" @@ -735,7 +735,11 @@ def plot_standard_performance_plots(args): "fakeorduplrate": ["pt", "ptlow", "ptmtv", "eta", "phi"], "avgOTlen": ["eta"], } - if (args.jet_branches): variables["eff"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "jet_eta", "jet_phi", "jet_pt"] + if (args.jet_branches): + variables["eff"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "genJetEta", "genJetPhi", "genJetPt"] + variables["duplrate"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "deltaR"] + variables["fakerate"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "deltaR"] + variables["fakeorduplrate"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "deltaR"] sels = { "eff": ["base", "loweta"], "fakerate": ["none"], @@ -757,9 +761,9 @@ def plot_standard_performance_plots(args): xcoarses["deltaEta"] = [False, True] xcoarses["deltaPhi"] = [False, True] xcoarses["deltaR"] = [False, True] - xcoarses["jet_eta"] = [False, True] - xcoarses["jet_phi"] = [False, True] - xcoarses["jet_pt"] = [False, True] + xcoarses["genJetEta"] = [False, True] + xcoarses["genJetPhi"] = [False, True] + xcoarses["genJetPt"] = [False, True] types = objecttype_choices breakdowns = { diff --git a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc index af1ab9d3f9edc..0dbf224bcc8f7 100644 --- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc +++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc @@ -725,52 +725,20 @@ void bookEfficiencySet(SimTrackSetDefinition& effset) { return ana.tx.getBranchLazy>(category_name + "_ie_numer_deltaR"); }); - // Lines for jet_eta - ana.tx.createBranch>(category_name + "_ef_denom_jet_eta"); - ana.histograms.addVecHistogram(category_name + "_ef_denom_jet_eta", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_denom_jet_eta"); + // Lines for genJetPt + ana.tx.createBranch>(category_name + "_ef_denom_genJetPt"); + ana.histograms.addVecHistogram(category_name + "_ef_denom_genJetPt", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_denom_genJetPt"); }); - ana.tx.createBranch>(category_name + "_ef_numer_jet_eta"); - ana.histograms.addVecHistogram(category_name + "_ef_numer_jet_eta", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_numer_jet_eta"); + ana.tx.createBranch>(category_name + "_ef_numer_genJetPt"); + ana.histograms.addVecHistogram(category_name + "_ef_numer_genJetPt", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_numer_genJetPt"); }); - ana.tx.createBranch>(category_name + "_ie_numer_jet_eta"); - ana.histograms.addVecHistogram(category_name + "_ie_numer_jet_eta", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ie_numer_jet_eta"); - }); - - // Lines for jet_phi - ana.tx.createBranch>(category_name + "_ef_denom_jet_phi"); - ana.histograms.addVecHistogram(category_name + "_ef_denom_jet_phi", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_denom_jet_phi"); - }); - - ana.tx.createBranch>(category_name + "_ef_numer_jet_phi"); - ana.histograms.addVecHistogram(category_name + "_ef_numer_jet_phi", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_numer_jet_phi"); - }); - - ana.tx.createBranch>(category_name + "_ie_numer_jet_phi"); - ana.histograms.addVecHistogram(category_name + "_ie_numer_jet_phi", 180, -4.5, 4.5, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ie_numer_jet_phi"); - }); - - // Lines for jet_pt - ana.tx.createBranch>(category_name + "_ef_denom_jet_pt"); - ana.histograms.addVecHistogram(category_name + "_ef_denom_jet_pt", 50, 50, 1000, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_denom_jet_pt"); - }); - - ana.tx.createBranch>(category_name + "_ef_numer_jet_pt"); - ana.histograms.addVecHistogram(category_name + "_ef_numer_jet_pt", 50, 50, 1000, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ef_numer_jet_pt"); - }); - - ana.tx.createBranch>(category_name + "_ie_numer_jet_pt"); - ana.histograms.addVecHistogram(category_name + "_ie_numer_jet_pt", 50, 50, 1000, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_ie_numer_jet_pt"); + ana.tx.createBranch>(category_name + "_ie_numer_genJetPt"); + ana.histograms.addVecHistogram(category_name + "_ie_numer_genJetPt", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ie_numer_genJetPt"); }); // Moving the standard pT code up here for convenience @@ -901,11 +869,13 @@ void bookDuplicateRateSet(RecoTrackSetDefinition& DRset) { ana.tx.createBranch>(category_name + "_dr_denom_pt"); ana.tx.createBranch>(category_name + "_dr_denom_eta"); ana.tx.createBranch>(category_name + "_dr_denom_phi"); + ana.tx.createBranch>(category_name + "_dr_denom_deltaR"); // Numerator tracks' quantities ana.tx.createBranch>(category_name + "_dr_numer_pt"); ana.tx.createBranch>(category_name + "_dr_numer_eta"); ana.tx.createBranch>(category_name + "_dr_numer_phi"); + ana.tx.createBranch>(category_name + "_dr_numer_deltaR"); // Histogram utility object that is used to define the histograms ana.histograms.addVecHistogram(category_name + "_dr_denom_pt", getPtBounds(0), [&, category_name]() { @@ -923,6 +893,9 @@ void bookDuplicateRateSet(RecoTrackSetDefinition& DRset) { ana.histograms.addVecHistogram(category_name + "_dr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_dr_denom_phi"); }); + ana.histograms.addVecHistogram(category_name + "_dr_denom_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_dr_denom_deltaR"); + }); ana.histograms.addVecHistogram(category_name + "_dr_numer_pt", getPtBounds(0), [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_dr_numer_pt"); }); @@ -938,6 +911,9 @@ void bookDuplicateRateSet(RecoTrackSetDefinition& DRset) { ana.histograms.addVecHistogram(category_name + "_dr_numer_phi", 180, -M_PI, M_PI, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_dr_numer_phi"); }); + ana.histograms.addVecHistogram(category_name + "_dr_numer_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_dr_numer_deltaR"); + }); } //__________________________________________________________________________________________________________________________________________________________________________ @@ -955,11 +931,13 @@ void bookFakeOrDuplicateRateSet(RecoTrackSetDefinition& FDRset) { ana.tx.createBranch>(category_name + "_fdr_denom_pt"); ana.tx.createBranch>(category_name + "_fdr_denom_eta"); ana.tx.createBranch>(category_name + "_fdr_denom_phi"); + ana.tx.createBranch>(category_name + "_fdr_denom_deltaR"); // Numerator tracks' quantities ana.tx.createBranch>(category_name + "_fdr_numer_pt"); ana.tx.createBranch>(category_name + "_fdr_numer_eta"); ana.tx.createBranch>(category_name + "_fdr_numer_phi"); + ana.tx.createBranch>(category_name + "_fdr_numer_deltaR"); // Histogram utility object that is used to define the histograms ana.histograms.addVecHistogram(category_name + "_fdr_denom_pt", getPtBounds(0), [&, category_name]() { @@ -977,6 +955,9 @@ void bookFakeOrDuplicateRateSet(RecoTrackSetDefinition& FDRset) { ana.histograms.addVecHistogram(category_name + "_fdr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fdr_denom_phi"); }); + ana.histograms.addVecHistogram(category_name + "_fdr_denom_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_fdr_denom_deltaR"); + }); ana.histograms.addVecHistogram(category_name + "_fdr_numer_pt", getPtBounds(0), [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fdr_numer_pt"); }); @@ -992,6 +973,9 @@ void bookFakeOrDuplicateRateSet(RecoTrackSetDefinition& FDRset) { ana.histograms.addVecHistogram(category_name + "_fdr_numer_phi", 180, -M_PI, M_PI, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fdr_numer_phi"); }); + ana.histograms.addVecHistogram(category_name + "_fdr_numer_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_fdr_numer_deltaR"); + }); } //__________________________________________________________________________________________________________________________________________________________________________ @@ -1009,11 +993,13 @@ void bookFakeRateSet(RecoTrackSetDefinition& FRset) { ana.tx.createBranch>(category_name + "_fr_denom_pt"); ana.tx.createBranch>(category_name + "_fr_denom_eta"); ana.tx.createBranch>(category_name + "_fr_denom_phi"); + ana.tx.createBranch>(category_name + "_fr_denom_deltaR"); // Numerator tracks' quantities ana.tx.createBranch>(category_name + "_fr_numer_pt"); ana.tx.createBranch>(category_name + "_fr_numer_eta"); ana.tx.createBranch>(category_name + "_fr_numer_phi"); + ana.tx.createBranch>(category_name + "_fr_numer_deltaR"); // Histogram utility object that is used to define the histograms ana.histograms.addVecHistogram(category_name + "_fr_denom_pt", getPtBounds(0), [&, category_name]() { @@ -1028,6 +1014,12 @@ void bookFakeRateSet(RecoTrackSetDefinition& FRset) { ana.histograms.addVecHistogram(category_name + "_fr_denom_eta", 180, -4.5, 4.5, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fr_denom_eta"); }); + ana.histograms.addVecHistogram(category_name + "_fr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_fr_denom_phi"); + }); + ana.histograms.addVecHistogram(category_name + "_fr_denom_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_fr_denom_deltaR"); + }); ana.histograms.addVecHistogram(category_name + "_fr_numer_phi", 180, -M_PI, M_PI, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fr_numer_phi"); }); @@ -1043,8 +1035,8 @@ void bookFakeRateSet(RecoTrackSetDefinition& FRset) { ana.histograms.addVecHistogram(category_name + "_fr_numer_eta", 180, -4.5, 4.5, [&, category_name]() { return ana.tx.getBranchLazy>(category_name + "_fr_numer_eta"); }); - ana.histograms.addVecHistogram(category_name + "_fr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() { - return ana.tx.getBranchLazy>(category_name + "_fr_denom_phi"); + ana.histograms.addVecHistogram(category_name + "_fr_numer_deltaR", 50, 0, 0.1, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_fr_numer_deltaR"); }); } @@ -1054,6 +1046,30 @@ void bookFakeRateSet(RecoTrackSetDefinition& FRset) { // ---------------------------------------------------------=============================================------------------------------------------------------------------- // ---------------------------------------------------------=============================================------------------------------------------------------------------- +float finddRTemp(float eta, + float phi, + std::vector genJetPt, + std::vector genJetEta, + std::vector genJetPhi) { + float dEtaj = 999; + float dPhij = 999; + float dRj2 = 999; + float dRTemp2 = 999; + float genJetPtj; + for (unsigned int ijet = 0; ijet < genJetPhi.size(); ++ijet) { + genJetPtj = genJetPt[ijet]; + if (genJetPtj > 1000) { + dEtaj = eta - genJetEta[ijet]; + dPhij = std::acos(std::cos(phi - genJetPhi[ijet])); + dRj2 = std::pow(dEtaj,2) + std::pow(dPhij,2); + if (dRj2 < dRTemp2) { + dRTemp2 = dRj2; + } + } + } + return std::sqrt(dRTemp2); +} + //__________________________________________________________________________________________________________________________________________________________________________ void fillEfficiencySets(std::vector& effsets) { std::vector const& pt = lstEff.getVF("sim_pt"); @@ -1068,33 +1084,34 @@ void fillEfficiencySets(std::vector& effsets) { std::vector const& vtx_z = lstEff.getVF("sim_vz"); if (ana.jet_branches) { - std::vector const& deltaEta = lstEff.getVF("sim_deltaEta"); - std::vector const& deltaPhi = lstEff.getVF("sim_deltaPhi"); - std::vector const& deltaR = lstEff.getVF("sim_deltaR"); - std::vector const& jet_eta = lstEff.getVF("sim_jet_eta"); - std::vector const& jet_phi = lstEff.getVF("sim_jet_phi"); - std::vector const& jet_pt = lstEff.getVF("sim_jet_pt"); - + std::vector const& deltaEta = lstEff.getVF("sim_genjet_deltaEta"); + std::vector const& deltaPhi = lstEff.getVF("sim_genjet_deltaPhi"); + std::vector const& deltaR = lstEff.getVF("sim_genjet_deltaR"); + std::vector const& genJetIdx = lstEff.getVI("sim_genjet_idx"); + std::vector const& genJetPt = lstEff.getVF("genjet_pt"); + float genJetPt_isimtrk = 0; + for (auto& effset : effsets) { for (unsigned int isimtrk = 0; isimtrk < lstEff.getVF("sim_pt").size(); ++isimtrk) { - fillEfficiencySet(isimtrk, - effset, - pt.at(isimtrk), - eta.at(isimtrk), - dz.at(isimtrk), - dxy.at(isimtrk), - phi.at(isimtrk), - pdgidtrk.at(isimtrk), - q.at(isimtrk), - vtx_x.at(isimtrk), - vtx_y.at(isimtrk), - vtx_z.at(isimtrk), - deltaEta.at(isimtrk), - deltaPhi.at(isimtrk), - deltaR.at(isimtrk), - jet_eta.at(isimtrk), - jet_phi.at(isimtrk), - jet_pt.at(isimtrk)); + genJetPt_isimtrk = genJetPt.at(genJetIdx.at(isimtrk)); + if (genJetPt_isimtrk > 1000){ + fillEfficiencySet(isimtrk, + effset, + pt.at(isimtrk), + eta.at(isimtrk), + dz.at(isimtrk), + dxy.at(isimtrk), + phi.at(isimtrk), + pdgidtrk.at(isimtrk), + q.at(isimtrk), + vtx_x.at(isimtrk), + vtx_y.at(isimtrk), + vtx_z.at(isimtrk), + deltaEta.at(isimtrk), + deltaPhi.at(isimtrk), + deltaR.at(isimtrk), + genJetPt_isimtrk); + } } } } else { @@ -1133,9 +1150,7 @@ void fillEfficiencySet(int isimtrk, float deltaEta, float deltaPhi, float deltaR, - float jet_eta, - float jet_phi, - float jet_pt) { + float genJetPt) { //========================================================= // NOTE: The following is not applied as the LSTNtuple no longer writes this. // const int &bunch = lstEff.getVI("sim_bunchCrossing")[isimtrk]; @@ -1173,103 +1188,87 @@ void fillEfficiencySet(int isimtrk, const float vtx_z_thresh = 30; const float vtx_perp_thresh = 2.5; - if (pt > 0 && jet_eta < 140 && jet_eta > -140 && (jet_eta > -999 && deltaEta > -999)) { - // N minus eta cut - if (pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh and std::abs(vtx_perp) < vtx_perp_thresh) { - // vs. eta plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_eta", eta); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_eta", eta); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_eta", eta); - } + // N minus eta cut + if (pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh and std::abs(vtx_perp) < vtx_perp_thresh) { + // vs. eta plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_eta", eta); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_eta", eta); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_eta", eta); + } - // N minus pt cut - if (std::abs(eta) < ana.eta_cut and std::abs(vtx_z) < vtx_z_thresh and std::abs(vtx_perp) < vtx_perp_thresh) { - // vs. pt plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_pt", pt); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_pt", pt); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_pt", pt); - } + // N minus pt cut + if (std::abs(eta) < ana.eta_cut and std::abs(vtx_z) < vtx_z_thresh and std::abs(vtx_perp) < vtx_perp_thresh) { + // vs. pt plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_pt", pt); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_pt", pt); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_pt", pt); + } - // N minus dxy cut - if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh) { - // vs. dxy plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_dxy", dxy); - ana.tx.pushbackToBranch(category_name + "_ef_denom_vxy", vtx_perp); - if (pass) { - ana.tx.pushbackToBranch(category_name + "_ef_numer_dxy", dxy); - ana.tx.pushbackToBranch(category_name + "_ef_numer_vxy", vtx_perp); - } else { - ana.tx.pushbackToBranch(category_name + "_ie_numer_dxy", dxy); - ana.tx.pushbackToBranch(category_name + "_ie_numer_vxy", vtx_perp); - } + // N minus dxy cut + if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh) { + // vs. dxy plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_dxy", dxy); + ana.tx.pushbackToBranch(category_name + "_ef_denom_vxy", vtx_perp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_ef_numer_dxy", dxy); + ana.tx.pushbackToBranch(category_name + "_ef_numer_vxy", vtx_perp); + } else { + ana.tx.pushbackToBranch(category_name + "_ie_numer_dxy", dxy); + ana.tx.pushbackToBranch(category_name + "_ie_numer_vxy", vtx_perp); } + } - // N minus dz cut - if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_perp) < vtx_perp_thresh) { - // vs. dz plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_dz", dz); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_dz", dz); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_dz", dz); - } + // N minus dz cut + if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_perp) < vtx_perp_thresh) { + // vs. dz plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_dz", dz); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_dz", dz); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_dz", dz); + } - // All phase-space cuts - if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh and - std::abs(vtx_perp) < vtx_perp_thresh) { - // vs. Phi plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_phi", phi); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_phi", phi); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_phi", phi); - - // vs. deltaEta plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaEta", deltaEta); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaEta", deltaEta); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaEta", deltaEta); - - // vs. deltaPhi plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaPhi", deltaPhi); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaPhi", deltaPhi); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaPhi", deltaPhi); - - // vs. deltaR plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaR", deltaR); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaR", deltaR); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaR", deltaR); - - // vs. jet_eta plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_jet_eta", jet_eta); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_jet_eta", jet_eta); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_jet_eta", jet_eta); - - // vs. jet_phi plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_jet_phi", jet_phi); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_jet_phi", jet_phi); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_jet_phi", jet_phi); - - // vs. jet_pt plot - ana.tx.pushbackToBranch(category_name + "_ef_denom_jet_pt", jet_pt); - if (pass) - ana.tx.pushbackToBranch(category_name + "_ef_numer_jet_pt", jet_pt); - else - ana.tx.pushbackToBranch(category_name + "_ie_numer_jet_pt", jet_pt); - } + // All phase-space cuts + if (std::abs(eta) < ana.eta_cut and pt > ana.pt_cut and std::abs(vtx_z) < vtx_z_thresh and + std::abs(vtx_perp) < vtx_perp_thresh) { + // vs. Phi plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_phi", phi); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_phi", phi); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_phi", phi); + + // vs. deltaEta plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaEta", deltaEta); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaEta", deltaEta); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaEta", deltaEta); + + // vs. deltaPhi plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaPhi", deltaPhi); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaPhi", deltaPhi); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaPhi", deltaPhi); + + // vs. deltaR plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_deltaR", deltaR); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_deltaR", deltaR); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_deltaR", deltaR); + + // vs. genJetPt plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_genJetPt", genJetPt); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_genJetPt", genJetPt); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_genJetPt", genJetPt); } } @@ -1414,6 +1413,18 @@ void fillFakeRateSet(int itc, RecoTrackSetDefinition& FRset, float pt, float eta if (pass) ana.tx.pushbackToBranch(category_name + "_fr_numer_phi", phi); } + if (ana.jet_branches) { + std::vector genJetPt = lstEff.getVF("genjet_pt"); + std::vector genJetEta = lstEff.getVF("genjet_eta"); + std::vector genJetPhi = lstEff.getVF("genjet_phi"); + + float dRTemp = finddRTemp(eta, phi, genJetPt, genJetEta, genJetPhi); + + ana.tx.pushbackToBranch(category_name + "_fr_denom_deltaR", dRTemp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_fr_numer_deltaR", dRTemp); + } + } } //__________________________________________________________________________________________________________________________________________________________________________ @@ -1452,6 +1463,18 @@ void fillDuplicateRateSet(int itc, RecoTrackSetDefinition& DRset, float pt, floa if (pass) ana.tx.pushbackToBranch(category_name + "_dr_numer_phi", phi); } + if (ana.jet_branches) { + std::vector genJetPt = lstEff.getVF("genjet_pt"); + std::vector genJetEta = lstEff.getVF("genjet_eta"); + std::vector genJetPhi = lstEff.getVF("genjet_phi"); + + float dRTemp = finddRTemp(eta, phi, genJetPt, genJetEta, genJetPhi); + + ana.tx.pushbackToBranch(category_name + "_dr_denom_deltaR", dRTemp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_dr_numer_deltaR", dRTemp); + } + } } //__________________________________________________________________________________________________________________________________________________________________________ @@ -1490,4 +1513,16 @@ void fillFakeOrDuplicateRateSet(int itc, RecoTrackSetDefinition& FDRset, float p if (pass) ana.tx.pushbackToBranch(category_name + "_fdr_numer_phi", phi); } + if (ana.jet_branches) { + std::vector genJetPt = lstEff.getVF("genjet_pt"); + std::vector genJetEta = lstEff.getVF("genjet_eta"); + std::vector genJetPhi = lstEff.getVF("genjet_phi"); + + float dRTemp = finddRTemp(eta, phi, genJetPt, genJetEta, genJetPhi); + + ana.tx.pushbackToBranch(category_name + "_fdr_denom_deltaR", dRTemp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_fdr_numer_deltaR", dRTemp); + } + } } diff --git a/RecoTracker/LSTCore/standalone/efficiency/src/performance.h b/RecoTracker/LSTCore/standalone/efficiency/src/performance.h index 97e64cb751456..b16f1b75ba154 100644 --- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.h +++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.h @@ -34,9 +34,7 @@ void fillEfficiencySet(int isimtrk, float deltaEta, float deltaPhi, float deltaR, - float jet_eta, - float jet_phi, - float jet_pt); + float genJetPt); void fillEfficiencySet(int isimtrk, SimTrackSetDefinition& effset, float pt, diff --git a/Validation/RecoTrack/plugins/TrackingNtuple.cc b/Validation/RecoTrack/plugins/TrackingNtuple.cc index 9440b9d47ad88..97a6b9e65ccc0 100644 --- a/Validation/RecoTrack/plugins/TrackingNtuple.cc +++ b/Validation/RecoTrack/plugins/TrackingNtuple.cc @@ -93,6 +93,12 @@ #include "RecoTracker/FinalTrackSelectors/interface/getBestVertex.h" +// GenJet headers +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" +#include "DataFormats/JetReco/interface/GenJet.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" + #include #include #include @@ -628,6 +634,9 @@ class TrackingNtuple : public edm::one::EDAnalyzer { void fillTrackingVertices(const TrackingVertexRefVector& trackingVertices, const TrackingParticleRefKeyToIndex& tpKeyToIndex); + void fillGenJets(const edm::Event& iEvent, + const TrackingParticleRefVector& tpCollection); + struct SimHitData { std::vector matchingSimHit; std::vector chargeFraction; @@ -705,6 +714,8 @@ class TrackingNtuple : public edm::one::EDAnalyzer { HistoryBase tracer_; ParametersDefinerForTP parametersDefiner_; + const edm::EDGetTokenT tok_jets_; + TTree* t; // DetId branches @@ -1415,6 +1426,18 @@ class TrackingNtuple : public edm::one::EDAnalyzer { std::vector> simvtx_sourceSimIdx; // second index runs through source TrackingParticles std::vector> simvtx_daughterSimIdx; // second index runs through daughter TrackingParticles std::vector simpv_idx; + std::vector sim_genjet_deltaEta; // distance to closest GenJet + std::vector sim_genjet_deltaPhi; + std::vector sim_genjet_deltaR; + std::vector sim_genjet_idx; // index of GenJet for each Sim track + + //////////////////// + // GenJets + std::vector genjet_pt; + std::vector genjet_eta; + std::vector genjet_phi; + std::vector genjet_invisible_energy; + std::vector genjet_auxiliary_energy; }; // @@ -1473,7 +1496,8 @@ TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig) keepEleSimHits_(iConfig.getUntrackedParameter("keepEleSimHits")), saveSimHitsP3_(iConfig.getUntrackedParameter("saveSimHitsP3")), simHitBySignificance_(iConfig.getUntrackedParameter("simHitBySignificance")), - parametersDefiner_(iConfig.getUntrackedParameter("beamSpot"), consumesCollector()) { + parametersDefiner_(iConfig.getUntrackedParameter("beamSpot"), consumesCollector()), + tok_jets_(consumes(iConfig.getParameter("jetSource"))) { if (includeSeeds_) { seedTokens_ = edm::vector_transform(iConfig.getUntrackedParameter>("seedTracks"), @@ -2050,9 +2074,20 @@ TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig) t->Branch("simvtx_z", &simvtx_z); t->Branch("simvtx_sourceSimIdx", &simvtx_sourceSimIdx); t->Branch("simvtx_daughterSimIdx", &simvtx_daughterSimIdx); + t->Branch("sim_genjet_deltaEta", &sim_genjet_deltaEta); + t->Branch("sim_genjet_deltaPhi", &sim_genjet_deltaPhi); + t->Branch("sim_genjet_deltaR", &sim_genjet_deltaR); t->Branch("simpv_idx", &simpv_idx); + // GenJets + t->Branch("genjet_pt", &genjet_pt); + t->Branch("genjet_eta", &genjet_eta); + t->Branch("genjet_phi", &genjet_phi); + t->Branch("genjet_invisible_energy", &genjet_invisible_energy); + t->Branch("genjet_auxiliary_energy", &genjet_auxiliary_energy); + t->Branch("sim_genjet_idx", &sim_genjet_idx); + //t->Branch("" , &); } @@ -2462,6 +2497,17 @@ void TrackingNtuple::clearVariables() { simvtx_sourceSimIdx.clear(); simvtx_daughterSimIdx.clear(); simpv_idx.clear(); + sim_genjet_deltaEta.clear(); + sim_genjet_deltaPhi.clear(); + sim_genjet_deltaR.clear(); + sim_genjet_idx.clear(); + + // GenJets + genjet_pt.clear(); + genjet_eta.clear(); + genjet_phi.clear(); + genjet_invisible_energy.clear(); + genjet_auxiliary_energy.clear(); } // ------------ method called for each event ------------ @@ -2735,9 +2781,53 @@ void TrackingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iS // tracking vertices fillTrackingVertices(tvRefs, tpKeyToIndex); + // GenJets + fillGenJets(iEvent, tpCollection); + t->Fill(); } +void TrackingNtuple::fillGenJets(const edm::Event& iEvent, + const TrackingParticleRefVector& tpCollection) { + auto const& genJets = iEvent.get(tok_jets_); + for (auto const& jet : genJets) { + genjet_pt.push_back(jet.pt()); + genjet_eta.push_back(jet.eta()); + genjet_phi.push_back(jet.phi()); + genjet_invisible_energy.push_back(jet.invisibleEnergy()); + genjet_auxiliary_energy.push_back(jet.auxiliaryEnergy()); + } + for (const TrackingParticleRef& tp : tpCollection) { + float sim_eta = tp->eta(); + float sim_phi = tp->phi(); + float dEtaj = 999; + float dPhij = 999; + float dRj2 = 999; + float dEtaTemp = 999; + float dPhiTemp = 999; + float dRTemp2 = 999; + int jet_idx = 999; + int counter = 0; + + for (auto const& jet : genJets) { + dEtaj = sim_eta - jet.eta(); + dPhij = reco::deltaPhi(sim_phi, jet.phi()); + dRj2 = std::pow(dEtaj,2) + std::pow(dPhij,2); + if(dRj2 < dRTemp2) { + dEtaTemp = dEtaj; + dPhiTemp = dPhij; + dRTemp2 = dRj2; + jet_idx = counter; + } + counter = counter + 1; + } + sim_genjet_deltaEta.push_back(dEtaTemp); + sim_genjet_deltaPhi.push_back(dPhiTemp); + sim_genjet_deltaR.push_back(std::sqrt(dRTemp2)); + sim_genjet_idx.push_back(jet_idx); + } +} + void TrackingNtuple::fillBeamSpot(const reco::BeamSpot& bs) { bsp_x = bs.x0(); bsp_y = bs.y0(); @@ -4721,6 +4811,7 @@ void TrackingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descriptio desc.addUntracked("keepEleSimHits", false); desc.addUntracked("saveSimHitsP3", false); desc.addUntracked("simHitBySignificance", false); + desc.add("jetSource", edm::InputTag("ak4GenJets")); descriptions.add("trackingNtuple", desc); }