diff --git a/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets.py b/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets.py index e3876fc6ecb33..b55f19d7620cf 100644 --- a/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets.py +++ b/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets.py @@ -13,12 +13,12 @@ import numpy as np # Load existing tree -file = TFile("/data2/segmentlinking/CMSSW_12_2_0_pre2/trackingNtuple_ttbar_PU200.root") -# file = TFile("trackingNtuple100.root") +# file = TFile("/data2/segmentlinking/CMSSW_12_2_0_pre2/trackingNtuple_ttbar_PU200.root") +file = TFile("trackingNtuple_100_GenJet.root") old_tree = file["trackingNtuple"]["tree"] # Create a new ROOT file to store the new TTree -new_file = ROOT.TFile("new_tree_temp.root", "RECREATE") +new_file = ROOT.TFile("new_tree_100_temp.root", "RECREATE") # Create a new subdirectory in the new file new_directory = new_file.mkdir("trackingNtuple") diff --git a/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets2.py b/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets2.py new file mode 100644 index 0000000000000..9ac4f747033f4 --- /dev/null +++ b/RecoTracker/LSTCore/standalone/analysis/jets/reformat_jets2.py @@ -0,0 +1,104 @@ +########################################################## +# +# Adds deltaR branches to the trackingNtuple.root file +# using GenJets. +# +########################################################## + +import matplotlib.pyplot as plt +import ROOT +from ROOT import TFile +import numpy as np + +# Load existing tree +# file = TFile("/data2/segmentlinking/CMSSW_12_2_0_pre2/trackingNtuple_ttbar_PU200.root") +file = TFile("trackingNtuple_100_GenJet.root") +old_tree = file["trackingNtuple"]["tree"] + +# Create a new ROOT file to store the new TTree +new_file = ROOT.TFile("new_tree_100_GenJet2.root", "RECREATE") + +# Create a new subdirectory in the new file +new_directory = new_file.mkdir("trackingNtuple") + +# Change the current directory to the new subdirectory +new_directory.cd() + +# Create a new TTree with the same structure as the old one but empty +new_tree = old_tree.CloneTree(0) + +# Account for bug in 12_2_X branch +new_tree.SetBranchStatus("ph2_bbxi", False) + +# Create a variable to hold the new leaves' data (a list of floats) +new_leaf_deltaEta = ROOT.std.vector('float')() +new_leaf_deltaPhi = ROOT.std.vector('float')() +new_leaf_deltaR = ROOT.std.vector('float')() + +# Create a new branch in the tree +new_tree.Branch("sim_deltaEta", new_leaf_deltaEta) +new_tree.Branch("sim_deltaPhi", new_leaf_deltaPhi) +new_tree.Branch("sim_deltaR", new_leaf_deltaR) + +# Loop over entries in the old tree +for ind in range(old_tree.GetEntries()): + old_tree.GetEntry(ind) + + # Clear the vector to start fresh for this entry + new_leaf_deltaEta.clear() + new_leaf_deltaPhi.clear() + new_leaf_deltaR.clear() + + # Creates the lists that will fill the leaves + simPt = old_tree.sim_pt + simEta = old_tree.sim_eta + simPhi = old_tree.sim_phi + simLen = len(simPt) + + print(old_tree.Print()) + + genJetsPt = old_tree.genJetPt + genJetsEta = old_tree.genJetEta + genJetsPhi = old_tree.genJetPhi + genJetLen = len(genJetsPt) + + # Declare arrays (all entries will be filled with non-dummy values) + deltaEtas = np.ones(simLen)*-999 + deltaPhis = np.ones(simLen)*-999 + deltaRs = np.ones(simLen)*-999 + + for i in range(len(simPt)): + dRTemp = 999 + dPhiTemp = 999 + dEtaTemp = 999 + for j in range(genJetLen): + dEtaj = simEta[i] - genJetsEta[j] + dPhij = np.arccos(np.cos(simPhi[i] - genJetsPhi[j])) + dRj = np.sqrt(dEtaj**2 + dPhij**2) + + # Selects smallest dR, corresponding to closest jet + if(dRj < dRTemp): + dRTemp = dRj + dPhiTemp = dPhij + dEtaTemp = dEtaj + deltaRs[i] = dRTemp + deltaPhis[i] = dPhiTemp + deltaEtas[i] = dEtaTemp + + + # Add the list elements to the vector + for value in deltaEtas: + new_leaf_deltaEta.push_back(value) + for value in deltaPhis: + new_leaf_deltaPhi.push_back(value) + for value in deltaRs: + new_leaf_deltaR.push_back(value) + + # Fill the tree with the new values + new_tree.Fill() + +# Write the tree back to the file +new_tree.Write() +new_file.Close() +file.Close() + diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc index 6428196e2c096..eda222295f059 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc @@ -195,9 +195,9 @@ 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>("genJetEta"); + ana.tx->createBranch>("genJetPhi"); + ana.tx->createBranch>("genJetPt"); } //________________________________________________________________________________________________________________________________ @@ -672,16 +672,16 @@ unsigned int setSimTrackContainerBranches(LSTEvent* event) { 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"); + auto const& trk_genJetEta = trk.getVF("genJetEta"); + auto const& trk_genJetPhi = trk.getVF("genJetPhi"); + auto const& trk_genJetPt = trk.getVF("genJetPt"); 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]); + ana.tx->pushbackToBranch("genJetEta", trk_genJetEta[isimtrk]); + ana.tx->pushbackToBranch("genJetPhi", trk_genJetPhi[isimtrk]); + ana.tx->pushbackToBranch("genJetPt", trk_genJetPt[isimtrk]); } // Fill the branch with simulated 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..0575b141d07a5 100644 --- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc +++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc @@ -725,52 +725,52 @@ 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 genJetEta + ana.tx.createBranch>(category_name + "_ef_denom_genJetEta"); + ana.histograms.addVecHistogram(category_name + "_ef_denom_genJetEta", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_denom_genJetEta"); }); - 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_genJetEta"); + ana.histograms.addVecHistogram(category_name + "_ef_numer_genJetEta", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_numer_genJetEta"); }); - 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"); + ana.tx.createBranch>(category_name + "_ie_numer_genJetEta"); + ana.histograms.addVecHistogram(category_name + "_ie_numer_genJetEta", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ie_numer_genJetEta"); }); - // 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"); + // Lines for genJetPhi + ana.tx.createBranch>(category_name + "_ef_denom_genJetPhi"); + ana.histograms.addVecHistogram(category_name + "_ef_denom_genJetPhi", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_denom_genJetPhi"); }); - 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 + "_ef_numer_genJetPhi"); + ana.histograms.addVecHistogram(category_name + "_ef_numer_genJetPhi", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_numer_genJetPhi"); }); - 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"); + ana.tx.createBranch>(category_name + "_ie_numer_genJetPhi"); + ana.histograms.addVecHistogram(category_name + "_ie_numer_genJetPhi", 180, -4.5, 4.5, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ie_numer_genJetPhi"); }); - // 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"); + // Lines for genJetPt + ana.tx.createBranch>(category_name + "_ef_denom_genJetPt"); + ana.histograms.addVecHistogram(category_name + "_ef_denom_genJetPt", 50, 50, 1000, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_denom_genJetPt"); }); - 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 + "_ef_numer_genJetPt"); + ana.histograms.addVecHistogram(category_name + "_ef_numer_genJetPt", 50, 50, 1000, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ef_numer_genJetPt"); }); - 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, 50, 1000, [&, category_name]() { + return ana.tx.getBranchLazy>(category_name + "_ie_numer_genJetPt"); }); // Moving the standard pT code up here for convenience @@ -901,11 +901,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 +925,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 +943,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 +963,14 @@ 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 +988,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 +1006,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 +1026,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 +1047,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 +1068,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"); }); } @@ -1071,30 +1096,32 @@ void fillEfficiencySets(std::vector& effsets) { 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& genJetEta = lstEff.getVF("genJetEta"); + std::vector const& genJetPhi = lstEff.getVF("genJetPhi"); + std::vector const& genJetPt = lstEff.getVF("genJetPt"); 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)); + if (genJetPt.at(isimtrk) > 100) { + 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), + genJetEta.at(isimtrk), + genJetPhi.at(isimtrk), + genJetPt.at(isimtrk)); + } } } } else { @@ -1133,9 +1160,9 @@ void fillEfficiencySet(int isimtrk, float deltaEta, float deltaPhi, float deltaR, - float jet_eta, - float jet_phi, - float jet_pt) { + float genJetEta, + float genJetPhi, + 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 +1200,101 @@ 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. genJetEta plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_genJetEta", genJetEta); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_genJetEta", genJetEta); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_genJetEta", genJetEta); + + // vs. genJetPhi plot + ana.tx.pushbackToBranch(category_name + "_ef_denom_genJetPhi", genJetPhi); + if (pass) + ana.tx.pushbackToBranch(category_name + "_ef_numer_genJetPhi", genJetPhi); + else + ana.tx.pushbackToBranch(category_name + "_ie_numer_genJetPhi", genJetPhi); + + // 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 +1439,33 @@ 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("genJetPt"); + std::vector genJetEta = lstEff.getVF("genJetEta"); + std::vector genJetPhi = lstEff.getVF("genJetPhi"); + + float dEtaj = 999; + float dPhij = 999; + float dRj = 999; + float dRTemp = 999; + float genJetPtj; + + for (unsigned int ijet; ijet < genJetPhi.size(); ++ijet) { + genJetPtj = genJetPt[ijet]; + if(genJetPtj > 100) { + dEtaj = eta - genJetEta[ijet]; + dPhij = std::acos(std::cos(phi - genJetPhi[ijet])); + dRj = std::sqrt(std::pow(dEtaj,2) + std::pow(dPhij,2)); + if(dRj < dRTemp) { + dRTemp = dRj; + } + } + } + ana.tx.pushbackToBranch(category_name + "_fr_denom_deltaR", dRTemp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_fr_numer_deltaR", dRTemp); + } + } } //__________________________________________________________________________________________________________________________________________________________________________ @@ -1452,6 +1504,34 @@ 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("genJetPt"); + std::vector genJetEta = lstEff.getVF("genJetEta"); + std::vector genJetPhi = lstEff.getVF("genJetPhi"); + + float dEtaj = 999; + float dPhij = 999; + float dRj = 999; + float dRTemp = 999; + float genJetPtj; + + for (unsigned int ijet; ijet < genJetPhi.size(); ++ijet) { + genJetPtj = genJetPt[ijet]; + if(genJetPtj > 100) { + dEtaj = eta - genJetEta[ijet]; + dPhij = std::acos(std::cos(phi - genJetPhi[ijet])); + dRj = std::sqrt(std::pow(dEtaj,2) + std::pow(dPhij,2)); + if(dRj < dRTemp) { + dRTemp = dRj; + } + } + } + ana.tx.pushbackToBranch(category_name + "_dr_denom_deltaR", dRTemp); + if (pass) { + ana.tx.pushbackToBranch(category_name + "_dr_numer_deltaR", dRTemp); + } + } } //__________________________________________________________________________________________________________________________________________________________________________ @@ -1490,4 +1570,31 @@ 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("genJetPt"); + std::vector genJetEta = lstEff.getVF("genJetEta"); + std::vector genJetPhi = lstEff.getVF("genJetPhi"); + + float dEtaj = 999; + float dPhij = 999; + float dRj = 999; + float dRTemp = 999; + float genJetPtj; + + for (unsigned int ijet; ijet < genJetPhi.size(); ++ijet) { + genJetPtj = genJetPt[ijet]; + if(genJetPtj > 100) { + dEtaj = eta - genJetEta[ijet]; + dPhij = std::acos(std::cos(phi - genJetPhi[ijet])); + dRj = std::sqrt(std::pow(dEtaj,2) + std::pow(dPhij,2)); + if(dRj < dRTemp) { + dRTemp = dRj; + } + } + } + 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..6092621bc1bc2 100644 --- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.h +++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.h @@ -34,9 +34,9 @@ void fillEfficiencySet(int isimtrk, float deltaEta, float deltaPhi, float deltaR, - float jet_eta, - float jet_phi, - float jet_pt); + float genJetEta, + float genJetPhi, + 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..4b69d355c8b6c 100644 --- a/Validation/RecoTrack/plugins/TrackingNtuple.cc +++ b/Validation/RecoTrack/plugins/TrackingNtuple.cc @@ -93,6 +93,13 @@ #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 @@ -627,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; @@ -705,6 +715,8 @@ class TrackingNtuple : public edm::one::EDAnalyzer { HistoryBase tracer_; ParametersDefinerForTP parametersDefiner_; + const edm::EDGetTokenT tok_jets_; + TTree* t; // DetId branches @@ -1295,7 +1307,7 @@ class TrackingNtuple : public edm::one::EDAnalyzer { ph2_radL; //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h std::vector ph2_bbxi; std::vector ph2_usedMask; - std::vector ph2_clustSize; + std::vector ph2_clustSize; //////////////////// // invalid (missing/inactive/etc) hits @@ -1415,6 +1427,17 @@ 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; + + //////////////////// + // GenJets + std::vector genjet_pt; + std::vector genjet_eta; + std::vector genjet_phi; + std::vector genjet_invisible_energy; + std::vector genjet_auxiliary_energy; + std::vector sim_deltaEta; + std::vector sim_deltaPhi; + std::vector sim_deltaR; }; // @@ -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"), @@ -1550,7 +1574,7 @@ TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig) index, consumes(mask.getUntrackedParameter("src"))); } } - + const bool tpRef = iConfig.getUntrackedParameter("trackingParticlesRef"); const auto tpTag = iConfig.getUntrackedParameter("trackingParticles"); if (tpRef) { @@ -2053,6 +2077,16 @@ TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig) 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_deltaEta", &sim_deltaEta); + t->Branch("sim_deltaPhi", &sim_deltaPhi); + t->Branch("sim_deltaR", &sim_deltaR); + //t->Branch("" , &); } @@ -2462,6 +2496,16 @@ void TrackingNtuple::clearVariables() { simvtx_sourceSimIdx.clear(); simvtx_daughterSimIdx.clear(); simpv_idx.clear(); + + // GenJets + genjet_pt.clear(); + genjet_eta.clear(); + genjet_phi.clear(); + genjet_invisible_energy.clear(); + genjet_auxiliary_energy.clear(); + sim_deltaEta.clear(); + sim_deltaPhi.clear(); + sim_deltaR.clear(); } // ------------ method called for each event ------------ @@ -2735,9 +2779,47 @@ 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 dRj = 999; + float dEtaTemp = 999; + float dPhiTemp = 999; + float dRTemp = 999; + for (auto const& jet : genJets) { + dEtaj = sim_eta - jet.eta(); + dPhij = std::acos(std::cos(sim_phi - jet.phi())); + dRj = std::sqrt(std::pow(dEtaj,2) + std::pow(dPhij,2)); + if(dRj < dRTemp) { + dEtaTemp = dEtaj; + dPhiTemp = dPhij; + dRTemp = dRj; + } + } + sim_deltaEta.push_back(dEtaTemp); + sim_deltaPhi.push_back(dPhiTemp); + sim_deltaR.push_back(dRTemp); + } +} + void TrackingNtuple::fillBeamSpot(const reco::BeamSpot& bs) { bsp_x = bs.x0(); bsp_y = bs.y0(); @@ -4721,8 +4803,9 @@ 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); } //define this as a plug-in -DEFINE_FWK_MODULE(TrackingNtuple); +DEFINE_FWK_MODULE(TrackingNtuple); \ No newline at end of file