diff --git a/analyzers/dataframe/FCCAnalyses/Analysis_FCChh.h b/analyzers/dataframe/FCCAnalyses/Analysis_FCChh.h index ca81379e71..0740c3fb2c 100644 --- a/analyzers/dataframe/FCCAnalyses/Analysis_FCChh.h +++ b/analyzers/dataframe/FCCAnalyses/Analysis_FCChh.h @@ -136,10 +136,12 @@ namespace AnalysisFCChh{ //btags ROOT::VecOps::RVec getJet_tag(ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec values, int algoIndex); ROOT::VecOps::RVec getBhadron(ROOT::VecOps::RVec truth_particles, ROOT::VecOps::RVec parent_ids); - ROOT::VecOps::RVec getChadron(ROOT::VecOps::RVec truth_particles, ROOT::VecOps::RVec parent_ids); - ROOT::VecOps::RVec get_tagged_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec values, int algoIndex); + ROOT::VecOps::RVec getChadron(ROOT::VecOps::RVec truth_particles, ROOT::VecOps::RVec parent_ids); + ROOT::VecOps::RVec get_tagged_jets(ROOT::VecOps::RVec reco_particles, ROOT::VecOps::RVec pids, ROOT::VecOps::RVec pids_rp_indices, ROOT::VecOps::RVec tag_values, int algoIndex); ROOT::VecOps::RVec get_untagged_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec values, int algoIndex); + + //tau jets ROOT::VecOps::RVec find_truth_matches(ROOT::VecOps::RVec truth_parts, ROOT::VecOps::RVec reco_particles, float dR_thres); ROOT::VecOps::RVec get_tau_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec tag_values, int algoIndex); diff --git a/analyzers/dataframe/src/Analysis_FCChh.cc b/analyzers/dataframe/src/Analysis_FCChh.cc index 67c8c8b5f6..b317b9920a 100644 --- a/analyzers/dataframe/src/Analysis_FCChh.cc +++ b/analyzers/dataframe/src/Analysis_FCChh.cc @@ -870,34 +870,62 @@ ROOT::VecOps::RVec AnalysisFCChh::getBhadron(ROOT::VecO } +//rewrite of functions to get tagged jets to work with updated edm4hep, https://github.com/key4hep/EDM4hep/pull/268 - -//return the full jets rather than the list of tags -ROOT::VecOps::RVec AnalysisFCChh::get_tagged_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec tag_values, int algoIndex){ +ROOT::VecOps::RVec AnalysisFCChh::get_tagged_jets(ROOT::VecOps::RVec reco_particles, + ROOT::VecOps::RVec pids, + ROOT::VecOps::RVec pids_rp_indices, + ROOT::VecOps::RVec tag_values, + int algoIndex){ - ROOT::VecOps::RVec tagged_jets; + ROOT::VecOps::RVec tagged_jets; - // std::cout << "running AnalysisFCChh::get_tagged_jets() on jet collection of size" << jets.size() << std::endl; + // make sure we have the right collections: every PID should have exactly one RP index + assert(pids.size() == pids_rp_indices.size()); + + + for (size_t pid_index = 0; pid_index < pids.size(); ++pid_index){ + + const auto tag = static_cast(tag_values[pids[pid_index].parameters_begin]); - for (size_t iJet = 0; iJet < jets.size(); ++iJet){ - // for (auto & jet : jets){ - std::cout << jet.particles_begin << " to " << jet.particles_end << std::endl; - // get the jet particle id index for the jet - const auto jetIDIndex = index[jets[iJet].particles_begin]; - // std::cout << "jet index = " << jetIDIndex << std::endl; - const auto jetID = pid[jetIDIndex]; - // get the tag value - const auto tag = static_cast(tag_values[jetID.parameters_begin]); // std::cout << "Tag = " << tag << std::endl; - // check if the tag satisfies what we want - if (tag & (1 << algoIndex)) { - tagged_jets.push_back(jets[iJet]); + + if (tag & (1 << algoIndex)) { + // std::cout << "Requested tag is true" << std::endl; + tagged_jets.push_back(reco_particles[pids_rp_indices[pid_index].index]); } - } + } return tagged_jets; } + +//return the full jets rather than the list of tags +// ROOT::VecOps::RVec AnalysisFCChh::get_tagged_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec tag_values, int algoIndex){ + +// ROOT::VecOps::RVec tagged_jets; + +// // std::cout << "running AnalysisFCChh::get_tagged_jets() on jet collection of size" << jets.size() << std::endl; + +// for (size_t iJet = 0; iJet < jets.size(); ++iJet){ +// // for (auto & jet : jets){ +// // std::cout << jet.particles_begin << " to " << jet.particles_end << std::endl; +// // get the jet particle id index for the jet +// const auto jetIDIndex = index[jets[iJet].particles_begin]; +// // std::cout << "jet index = " << jetIDIndex << std::endl; +// const auto jetID = pid[jetIDIndex]; +// // get the tag value +// const auto tag = static_cast(tag_values[jetID.parameters_begin]); +// // std::cout << "Tag = " << tag << std::endl; +// // check if the tag satisfies what we want +// if (tag & (1 << algoIndex)) { +// tagged_jets.push_back(jets[iJet]); +// } +// } + +// return tagged_jets; +// } + //same for tau tags: they are second entry in the ROOT::VecOps::RVec AnalysisFCChh::get_tau_jets(ROOT::VecOps::RVec jets, ROOT::VecOps::RVec index, ROOT::VecOps::RVec pid, ROOT::VecOps::RVec tag_values, int algoIndex){ diff --git a/examples/FCChh/ggHH_bbyy/analysis_stage1.py b/examples/FCChh/ggHH_bbyy/analysis_stage1.py index 6dc10346e6..845301345e 100644 --- a/examples/FCChh/ggHH_bbyy/analysis_stage1.py +++ b/examples/FCChh/ggHH_bbyy/analysis_stage1.py @@ -68,12 +68,11 @@ def analyzers(self, dframe): .Define("weight", "EventHeader.weight") ########################################### PHOTONS ########################################### - # .Alias("Photon0", "Photon#0.index") - # #.Alias("Photon0", "Photon_objIdx.index") + .Define("gamma", "FCCAnalyses::ReconstructedParticle::get(Photon_objIdx.index, ReconstructedParticles)") - .Define("sel_gamma", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(gamma)") - # .Define("sel_gamma_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_gamma)") - # .Define("sel_gamma", "AnalysisFCChh::SortParticleCollection(sel_gamma_unsort)") + .Define("selpt_gamma", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(gamma)") + .Define("sel_gamma_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_gamma)") + .Define("sel_gamma", "AnalysisFCChh::SortParticleCollection(sel_gamma_unsort)") .Define("ngamma", "FCCAnalyses::ReconstructedParticle::get_n(sel_gamma)") .Define("g1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_gamma)[0]") @@ -85,69 +84,58 @@ def analyzers(self, dframe): .Define("g2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_gamma)[1]") .Define("g2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_gamma)[1]") + #H(yy) if it exists, if there are no 2 selected photons, doesnt get filled + .Define("yy_pairs_unmerged", "AnalysisFCChh::getPairs(sel_gamma)") + .Define("yy_pairs", "AnalysisFCChh::merge_pairs(yy_pairs_unmerged)") + .Define("m_yy", "FCCAnalyses::ReconstructedParticle::get_mass(yy_pairs)") + ########################################### JETS ########################################### # jets after overlap removal is performed between jets and isolated electrons, muons and photons - #selected jets above a pT threshold of 30 GeV, eta < 4, tight ID - # .Define("selpt_jets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(Jet)") - # .Define("sel_jets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_jets)") - # .Define("sel_jets", "AnalysisFCChh::SortParticleCollection(sel_jets_unsort)") - # .Define("njets", "FCCAnalyses::ReconstructedParticle::get_n(sel_jets)") - # .Define("j1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[0]") - # .Define("j1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[0]") - # .Define("j1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[0]") - # .Define("j1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[0]") - - # .Define("j2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[1]") - # .Define("j2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[1]") - # .Define("j2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[1]") - # .Define("j2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[1]") - - - # # define an alias for muon index collection - # .Alias('Muon0', 'Muon_objIdx.index') - # # define the muon collection - # .Define( - # 'muons', - # 'ReconstructedParticle::get(Muon0, ReconstructedParticles)') - # # select muons on pT - # .Define('selected_muons', - # f'ReconstructedParticle::sel_pt({muon_pt})(muons)') - # # create branch with muon transverse momentum - # .Define('selected_muons_pt', - # 'ReconstructedParticle::get_pt(selected_muons)') - # # create branch with muon rapidity - # .Define('selected_muons_y', - # 'ReconstructedParticle::get_y(selected_muons)') - # # create branch with muon total momentum - # .Define('selected_muons_p', - # 'ReconstructedParticle::get_p(selected_muons)') - # # create branch with muon energy - # .Define('selected_muons_e', - # 'ReconstructedParticle::get_e(selected_muons)') - # # find zed candidates from di-muon resonances - # .Define( - # 'zed_leptonic', - # 'ReconstructedParticle::resonanceBuilder(91)(selected_muons)') - # # create branch with zed mass - # .Define('zed_leptonic_m', - # 'ReconstructedParticle::get_mass(zed_leptonic)') - # # create branch with zed transverse momenta - # .Define('zed_leptonic_pt', - # 'ReconstructedParticle::get_pt(zed_leptonic)') - # # calculate recoil of zed_leptonic - # .Define('zed_leptonic_recoil', - # 'ReconstructedParticle::recoilBuilder(240)(zed_leptonic)') - # # create branch with recoil mass - # .Define('zed_leptonic_recoil_m', - # 'ReconstructedParticle::get_mass(zed_leptonic_recoil)') - # # create branch with leptonic charge - # .Define('zed_leptonic_charge', - # 'ReconstructedParticle::get_charge(zed_leptonic)') - # # Filter at least one candidate - # .Filter('zed_leptonic_recoil_m.size()>0') + # selected jets above a pT threshold of 30 GeV, eta < 4, tight ID + .Define("selpt_jets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(Jet)") + .Define("sel_jets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_jets)") + .Define("sel_jets", "AnalysisFCChh::SortParticleCollection(sel_jets_unsort)") + .Define("njets", "FCCAnalyses::ReconstructedParticle::get_n(sel_jets)") + .Define("j1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[0]") + .Define("j1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[0]") + .Define("j1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[0]") + .Define("j1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[0]") + + .Define("j2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[1]") + .Define("j2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[1]") + .Define("j2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[1]") + .Define("j2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[1]") + + #b-tagged jets: + # .Alias("Jet3","Jet#3.index") + .Alias("Jet3","_Jet_particles.index") + #LOOSE WP : bit 1 = medium WP, bit 2 = tight WP + #b tagged jets + .Define("bjets", "AnalysisFCChh::get_tagged_jets(ReconstructedParticles, ParticleIDs, _ParticleIDs_particle, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl + # .Define("bjets", "AnalysisFCChh::get_tagged_jets(Jet, Jet3, ParticleIDs, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl + .Define("selpt_bjets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(bjets)") + .Define("sel_bjets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_bjets)") + .Define("sel_bjets", "AnalysisFCChh::SortParticleCollection(sel_bjets_unsort)") + .Define("nbjets", "FCCAnalyses::ReconstructedParticle::get_n(sel_bjets)") + .Define("b1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_bjets)[0]") + .Define("b1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_bjets)[0]") + .Define("b1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_bjets)[0]") + .Define("b1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_bjets)[0]") + .Define("b2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_bjets)[1]") + .Define("b2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_bjets)[1]") + .Define("b2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_bjets)[1]") + .Define("b2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_bjets)[1]") + + + # #H(bb) system - using the loose WP b-jets -> if the events has < 2 bjets these variables do not get filled! + # .Define("bb_pairs_unmerged", "AnalysisFCChh::getPairs(sel_bjets)") #currently gets only leading pT pair, as a RecoParticlePair + # #then merge the bb pair into one object and get its kinematic properties + # .Define("bb_pairs", "AnalysisFCChh::merge_pairs(bb_pairs_unmerged)") #merge into one object to access inv masses etc + # .Define("m_bb", "FCCAnalyses::ReconstructedParticle::get_mass(bb_pairs)") + ) return dframe2 @@ -158,7 +146,10 @@ def output(self): Output variables which will be saved to output root file. ''' branch_list = [ + 'nbjets', + # 'm_bb', 'ngamma', + 'm_yy', # 'selected_muons_pt', # 'selected_muons_y', # 'selected_muons_p',