diff --git a/RecoTracker/LSTCore/standalone/bin/sdl.cc b/RecoTracker/LSTCore/standalone/bin/sdl.cc index bdba8f9d5e44f..dd3acb041a135 100644 --- a/RecoTracker/LSTCore/standalone/bin/sdl.cc +++ b/RecoTracker/LSTCore/standalone/bin/sdl.cc @@ -51,6 +51,7 @@ int main(int argc, char **argv) { "o,output", "Output file name", cxxopts::value())( "N,nmatch", "N match for MTV-like matching", cxxopts::value()->default_value("9"))( "p,ptCut", "Min pT cut In GeV", cxxopts::value()->default_value("0.8"))( + "f,fmatch", "hit matching fraction for RECO-SIM matching", cxxopts::value()->default_value("0.75"))( "n,nevents", "N events to loop over", cxxopts::value()->default_value("-1"))( "x,event_index", "specific event index to process", cxxopts::value()->default_value("-1"))( "g,pdg_id", "The simhit pdgId match option", cxxopts::value()->default_value("0"))( @@ -60,7 +61,8 @@ int main(int argc, char **argv) { "w,write_ntuple", "Write Ntuple", cxxopts::value()->default_value("1"))( "s,streams", "Set number of streams", cxxopts::value()->default_value("1"))( "d,debug", "Run debug job. i.e. overrides output option to 'debug.root' and 'recreate's the file.")( - "l,lower_level", "write lower level objects ntuple results")("G,gnn_ntuple", "write gnn input variable ntuple")( + "O,optional_output", "Write optional objects in output LST ntuple")( + "l,lower_level", "Write lower level objects in output LST ntuple")( "j,nsplit_jobs", "Enable splitting jobs by N blocks (--job_index must be set)", cxxopts::value())( "I,job_index", "job_index of split jobs (--nsplit_jobs must be set. index starts from 0. i.e. 0, 1, 2, 3, etc...)", @@ -151,6 +153,10 @@ int main(int argc, char **argv) { // --nmatch ana.nmatch_threshold = result["nmatch"].as(); + //_______________________________________________________________________________ + // --fmatch + ana.fmatch_threshold = result["fmatch"].as(); + //_______________________________________________________________________________ // --nevents ana.n_events = result["nevents"].as(); @@ -227,25 +233,30 @@ int main(int argc, char **argv) { // --optimization //_______________________________________________________________________________ - // --lower_level - if (result.count("lower_level")) { - ana.do_lower_level = true; + // --optional_output + if (result.count("optional_output")) { + ana.do_optional_output = true; + if (not ana.do_write_ntuple) { + std::cout << options.help() << std::endl; + std::cout << "ERROR: option string --write_ntuple 1 and --optional_output must be set at the same time!" + << std::endl; + exit(1); + } } else { - ana.do_lower_level = false; + ana.do_optional_output = false; } //_______________________________________________________________________________ - // --gnn_ntuple - if (result.count("gnn_ntuple")) { - ana.gnn_ntuple = true; - // If one is not provided then throw error + // --lower_level + if (result.count("lower_level")) { + ana.do_lower_level = true; if (not ana.do_write_ntuple) { std::cout << options.help() << std::endl; - std::cout << "ERROR: option string --write_ntuple 1 and --gnn_ntuple must be set at the same time!" << std::endl; + std::cout << "ERROR: option string --write_ntuple 1 and --lower_level must be set at the same time!" << std::endl; exit(1); } } else { - ana.gnn_ntuple = false; + ana.do_lower_level = false; } //_______________________________________________________________________________ @@ -273,6 +284,7 @@ int main(int argc, char **argv) { std::cout << " ana.verbose: " << ana.verbose << std::endl; std::cout << " ana.ptCut: " << ana.ptCut << std::endl; std::cout << " ana.nmatch_threshold: " << ana.nmatch_threshold << std::endl; + std::cout << " ana.fmatch_threshold: " << ana.fmatch_threshold << std::endl; std::cout << " ana.tc_pls_triplets: " << ana.tc_pls_triplets << std::endl; std::cout << " ana.no_pls_dupclean: " << ana.no_pls_dupclean << std::endl; std::cout << "=========================================================" << std::endl; @@ -319,8 +331,8 @@ void run_sdl() { if (ana.do_write_ntuple) { createOutputBranches(); - if (ana.gnn_ntuple) { - createGnnNtupleBranches(); + if (ana.do_lower_level) { + createLowLevelBranches(); } } diff --git a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc index dfae0e4925024..3828e8554cbbf 100644 --- a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc +++ b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc @@ -139,6 +139,33 @@ std::vector getHitsFromT3(SDL::Event* event, unsigned in return {hits_0[0], hits_0[1], hits_1[0], hits_1[1], hits_2[0], hits_2[1]}; } +//____________________________________________________________________________________________ +std::vector getHitIdxsFromT3(SDL::Event* event, unsigned int T3) { + SDL::hitsBuffer& hitsInGPU = *(event->getHits()); + std::vector hits = getHitsFromT3(event, T3); + std::vector hitidxs; + for (auto& hit : hits) + hitidxs.push_back(hitsInGPU.idxs[hit]); + return hitidxs; +} + +//____________________________________________________________________________________________ +std::vector getModuleIdxsFromT3(SDL::Event* event, unsigned int T3) { + std::vector hits = getHitsFromT3(event, T3); + std::vector module_idxs; + SDL::hitsBuffer& hitsInGPU = *(event->getHits()); + for (auto& hitIdx : hits) { + module_idxs.push_back(hitsInGPU.moduleIndices[hitIdx]); + } + return module_idxs; +} + +//____________________________________________________________________________________________ +std::vector getHitTypesFromT3(SDL::Event* event, unsigned int T3) { + return {4, 4, 4, 4, 4, 4}; + ; +} + //____________________________________________________________________________________________ std::tuple, std::vector> getHitIdxsAndHitTypesFromT3( SDL::Event* event, unsigned T3) { diff --git a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h index 589f1e2d59cd2..0c5986a81546a 100644 --- a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h +++ b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h @@ -33,6 +33,9 @@ std::tuple, std::vector> getHitIdxsAndHi std::vector getLSsFromT3(SDL::Event* event, unsigned int T3); std::vector getMDsFromT3(SDL::Event* event, unsigned int T3); std::vector getHitsFromT3(SDL::Event* event, unsigned int T3); +std::vector getHitIdxsFromT3(SDL::Event* event, unsigned int T3); +std::vector getHitTypesFromT3(SDL::Event* event, unsigned int T3); +std::vector getModuleIdxsFromT3(SDL::Event* event, unsigned int T3); std::tuple, std::vector> getHitIdxsAndHitTypesFromT3( SDL::Event* event, unsigned T3); diff --git a/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h b/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h index 50ef9a1ff5a71..8d99dbe6c7797 100644 --- a/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h +++ b/RecoTracker/LSTCore/standalone/code/core/AnalysisConfig.h @@ -55,6 +55,9 @@ class AnalysisConfig { // nmatch threshold of the hits to compute matching for MTV-like plots int nmatch_threshold; + // min matching hit fraction to determine RECO-SIM match + float fmatch_threshold; + // verbose of the particles to compute efficincies on int verbose; @@ -111,12 +114,12 @@ class AnalysisConfig { // Boolean to trigger whether to write ntuple or not bool do_write_ntuple; + // Boolean to write optional output + bool do_optional_output; + // Boolean to write lower level objects bool do_lower_level; - // Boolean to write gnn ntuple - bool gnn_ntuple; - // String to hold the MAKETARGET setting from compile std::string compilation_target; diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc index 1a8593d12d113..8741252cdc297 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc @@ -283,7 +283,10 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hi //___________________________________________________________________________________________________________________________________________________________________________________________ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, - bool verbose) { + bool verbose, + float matchFraction, + bool writeMatchFraction, + std::vector &matched_fractions) { if (hitidxs.size() != hittypes.size()) { std::cout << "Error: matched_sim_trk_idxs() hitidxs and hittypes have different lengths" << std::endl; std::cout << "hitidxs.size(): " << hitidxs.size() << std::endl; @@ -414,6 +417,7 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, } int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits std::vector matched_sim_trk_idxs; + std::vector> matched_sim_idx_frac; for (auto &trkidx_perm : allperms) { std::vector counts; for (auto &unique_idx : unique_idxs) { @@ -425,15 +429,22 @@ std::vector matchedSimTrkIdxs(std::vector hitidxs, int trkidx = unique_idxs[rawidx]; if (trkidx < 0) continue; - if (counts[rawidx] > (((float)nhits_input) * 0.75)) - matched_sim_trk_idxs.push_back(trkidx); + if (counts[rawidx] > (((float)nhits_input) * matchFraction)) + matched_sim_idx_frac.push_back( + std::make_pair(trkidx, (nhits_input > 0 ? ((float)counts[rawidx]) / ((float)nhits_input) : -999.0))); maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end())); } - set s; - unsigned size = matched_sim_trk_idxs.size(); - for (unsigned i = 0; i < size; ++i) - s.insert(matched_sim_trk_idxs[i]); - matched_sim_trk_idxs.assign(s.begin(), s.end()); + sort(matched_sim_idx_frac.begin(), matched_sim_idx_frac.end(), [](auto &a, auto &b) { return a.second > b.second; }); + for (auto it = std::make_move_iterator(matched_sim_idx_frac.begin()), + end = std::make_move_iterator(matched_sim_idx_frac.end()); + it != end; + ++it) { + if (std::find(matched_sim_trk_idxs.begin(), matched_sim_trk_idxs.end(), it->first) == matched_sim_trk_idxs.end()) { + matched_sim_trk_idxs.push_back(std::move(it->first)); + if (writeMatchFraction) + matched_fractions.push_back(std::move(it->second)); + } + } return matched_sim_trk_idxs; } diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.h b/RecoTracker/LSTCore/standalone/code/core/trkCore.h index a407f303cf1a6..f6992103cff69 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.h +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.h @@ -28,9 +28,13 @@ float runpT3(SDL::Event *event); // --------------------- ======================== --------------------- +static std::vector empty_vector_float = std::vector(); std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, - bool verbose = false); + bool verbose = false, + float matchFraction = 0.75, + bool writeMatchFraction = false, + std::vector &matched_fractions = empty_vector_float); std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose = false); int getDenomSimTrkType(int isimtrk); int getDenomSimTrkType(std::vector simidxs); diff --git a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc index 35a7507770dcd..9f673156fca74 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc @@ -9,9 +9,12 @@ void createOutputBranches() { //________________________________________________________________________________________________________________________________ void fillOutputBranches(SDL::Event* event) { setOutputBranches(event); - setOptionalOutputBranches(event); - if (ana.gnn_ntuple) - setGnnNtupleBranches(event); + if (ana.do_optional_output) + setOptionalOutputBranches(event); + if (ana.do_lower_level) { + setLowLevelBranches(event); + setTripletOutputBranches(event); + } // Now actually fill the ttree ana.tx->fill(); @@ -50,7 +53,6 @@ void createRequiredOutputBranches() { //________________________________________________________________________________________________________________________________ void createOptionalOutputBranches() { -#ifdef CUT_VALUE_DEBUG // Event-wide branches // ana.tx->createBranch("evt_dummy"); @@ -74,7 +76,6 @@ void createOptionalOutputBranches() { ana.tx->createBranch>("pT5_score"); ana.tx->createBranch>("pT5_layer_binary"); ana.tx->createBranch>("pT5_moduleType_binary"); - ana.tx->createBranch>("pT5_matched_pt"); ana.tx->createBranch>("pT5_rzChiSquared"); ana.tx->createBranch>("pT5_rPhiChiSquared"); ana.tx->createBranch>("pT5_rPhiChiSquaredInwards"); @@ -87,12 +88,10 @@ void createOptionalOutputBranches() { ana.tx->createBranch>("pT3_eta"); ana.tx->createBranch>("pT3_phi"); ana.tx->createBranch>("pT3_score"); - ana.tx->createBranch>("pT3_foundDuplicate"); ana.tx->createBranch>>("pT3_matched_simIdx"); ana.tx->createBranch>>("pT3_hitIdxs"); ana.tx->createBranch>("pT3_pixelRadius"); ana.tx->createBranch>("pT3_pixelRadiusError"); - ana.tx->createBranch>>("pT3_matched_pt"); ana.tx->createBranch>("pT3_tripletRadius"); ana.tx->createBranch>("pT3_rPhiChiSquared"); ana.tx->createBranch>("pT3_rPhiChiSquaredInwards"); @@ -114,7 +113,6 @@ void createOptionalOutputBranches() { ana.tx->createBranch>("sim_T5_matched"); ana.tx->createBranch>("t5_isFake"); ana.tx->createBranch>("t5_isDuplicate"); - ana.tx->createBranch>("t5_foundDuplicate"); ana.tx->createBranch>("t5_pt"); ana.tx->createBranch>("t5_eta"); ana.tx->createBranch>("t5_phi"); @@ -123,8 +121,7 @@ void createOptionalOutputBranches() { ana.tx->createBranch>>("t5_matched_simIdx"); ana.tx->createBranch>("t5_moduleType_binary"); ana.tx->createBranch>("t5_layer_binary"); - ana.tx->createBranch>("t5_matched_pt"); - ana.tx->createBranch>("t5_partOfTC"); + ana.tx->createBranch>("t5_partOfPT5"); ana.tx->createBranch>("t5_innerRadius"); ana.tx->createBranch>("t5_outerRadius"); ana.tx->createBranch>("t5_bridgeRadius"); @@ -148,12 +145,10 @@ void createOptionalOutputBranches() { ana.tx->createBranch>("t5_occupancies"); ana.tx->createBranch("pT3_occupancies"); ana.tx->createBranch("pT5_occupancies"); - -#endif } //________________________________________________________________________________________________________________________________ -void createGnnNtupleBranches() { +void createLowLevelBranches() { // Mini Doublets ana.tx->createBranch>("MD_pt"); ana.tx->createBranch>("MD_eta"); @@ -173,28 +168,37 @@ void createGnnNtupleBranches() { ana.tx->createBranch>("MD_1_z"); // Line Segments + ana.tx->createBranch>("sim_LS_matched"); + ana.tx->createBranch>("LS_isFake"); + ana.tx->createBranch>("LS_isDuplicate"); + ana.tx->createBranch>("LS_nsimMatch"); + ana.tx->createBranch>>("LS_matched_frac"); + ana.tx->createBranch>>("LS_matched_simIdx"); ana.tx->createBranch>("LS_pt"); ana.tx->createBranch>("LS_eta"); ana.tx->createBranch>("LS_phi"); - ana.tx->createBranch>("LS_isFake"); ana.tx->createBranch>("LS_MD_idx0"); ana.tx->createBranch>("LS_MD_idx1"); - ana.tx->createBranch>("LS_sim_pt"); - ana.tx->createBranch>("LS_sim_eta"); - ana.tx->createBranch>("LS_sim_phi"); - ana.tx->createBranch>("LS_sim_pca_dxy"); - ana.tx->createBranch>("LS_sim_pca_dz"); - ana.tx->createBranch>("LS_sim_q"); - ana.tx->createBranch>("LS_sim_pdgId"); - ana.tx->createBranch>("LS_sim_event"); - ana.tx->createBranch>("LS_sim_bx"); - ana.tx->createBranch>("LS_sim_vx"); - ana.tx->createBranch>("LS_sim_vy"); - ana.tx->createBranch>("LS_sim_vz"); ana.tx->createBranch>("LS_isInTrueTC"); // TC's LS ana.tx->createBranch>>("tc_lsIdx"); + + // T3 branches + ana.tx->createBranch>("sim_T3_matched"); + ana.tx->createBranch>("t3_isFake"); + ana.tx->createBranch>("t3_isDuplicate"); + ana.tx->createBranch>("t3_pt"); + ana.tx->createBranch>("t3_eta"); + ana.tx->createBranch>("t3_phi"); + ana.tx->createBranch>("t3_circleRadius"); + ana.tx->createBranch>>("t3_hitIdxs"); + ana.tx->createBranch>>("t3_matched_simIdx"); + ana.tx->createBranch>("t3_moduleType_binary"); + ana.tx->createBranch>("t3_layer_binary"); + ana.tx->createBranch>("t3_partOfPT3"); + ana.tx->createBranch>("t3_partOfPT5"); + ana.tx->createBranch>("t3_partOfT5"); } //________________________________________________________________________________________________________________________________ @@ -257,15 +261,15 @@ void setOutputBranches(SDL::Event* event) { tc_matched_simIdx.push_back(simidx); // Loop over matched sim idx and increase counter of TC_matched - for (auto& idx : simidx) { + for (auto& sidx : simidx) { // NOTE Important to note that the idx of the std::vector<> is same // as the tracking-ntuple's sim track idx ONLY because event==0 and bunchCrossing==0 condition is applied!! // Also do not try to access beyond the event and bunchCrossing - if (idx < n_accepted_simtrk) { - sim_TC_matched.at(idx) += 1; - sim_TC_matched_mask.at(idx) |= (1 << type); + if (sidx < n_accepted_simtrk) { + sim_TC_matched.at(sidx) += 1; + sim_TC_matched_mask.at(sidx) |= (1 << type); } - sim_TC_matched_for_duplicate.at(idx) += 1; + sim_TC_matched_for_duplicate.at(sidx) += 1; } } @@ -294,14 +298,10 @@ void setOutputBranches(SDL::Event* event) { //________________________________________________________________________________________________________________________________ void setOptionalOutputBranches(SDL::Event* event) { -#ifdef CUT_VALUE_DEBUG - setPixelQuintupletOutputBranches(event); setQuintupletOutputBranches(event); setPixelTripletOutputBranches(event); setOccupancyBranches(event); - -#endif } //________________________________________________________________________________________________________________________________ @@ -397,7 +397,7 @@ void setPixelQuintupletOutputBranches(SDL::Event* event) { layer_binary |= (1 << (modulesInGPU.layers[module_idx[i]] + 6 * (modulesInGPU.subdets[module_idx[i]] == 4))); moduleType_binary |= (modulesInGPU.moduleType[module_idx[i]] << i); } - std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type); + std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type, ana.verbose, ana.fmatch_threshold); ana.tx->pushbackToBranch("pT5_isFake", static_cast(simidx.size() == 0)); ana.tx->pushbackToBranch("pT5_pt", pt); ana.tx->pushbackToBranch("pT5_eta", eta); @@ -485,6 +485,7 @@ void setQuintupletOutputBranches(SDL::Event* event) { ana.tx->pushbackToBranch("t5_rzChiSquared", quintupletsInGPU.rzChiSquared[quintupletIndex]); ana.tx->pushbackToBranch("t5_layer_binary", layer_binary); ana.tx->pushbackToBranch("t5_moduleType_binary", moduleType_binary); + ana.tx->pushbackToBranch("t5_partOfPT5", __H2F(quintupletsInGPU.partOfPT5[quintupletIndex])); t5_matched_simIdx.push_back(simidx); @@ -537,7 +538,7 @@ void setPixelTripletOutputBranches(SDL::Event* event) { std::vector hit_idx = getHitIdxsFrompT3(event, pT3); std::vector hit_type = getHitTypesFrompT3(event, pT3); - std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type); + std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type, ana.verbose, ana.fmatch_threshold); std::vector module_idx = getModuleIdxsFrompT3(event, pT3); int layer_binary = 1; int moduleType_binary = 0; @@ -580,7 +581,95 @@ void setPixelTripletOutputBranches(SDL::Event* event) { } //________________________________________________________________________________________________________________________________ -void setGnnNtupleBranches(SDL::Event* event) { +void setTripletOutputBranches(SDL::Event* event) { + SDL::tripletsBuffer& tripletsInGPU = *(event->getTriplets()); + SDL::objectRangesBuffer& rangesInGPU = (*event->getRanges()); + SDL::modulesBuffer& modulesInGPU = *(event->getModules()); + SDL::segmentsBuffer& segmentsInGPU = *(event->getSegments()); + SDL::hitsBuffer& hitsInGPU = *(event->getHits()); + std::vector> t3_matched_simIdx; + + for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < *(modulesInGPU.nLowerModules); ++lowerModuleIdx) { + int nTriplets = tripletsInGPU.nTriplets[lowerModuleIdx]; + for (unsigned int idx = 0; idx < nTriplets; idx++) { + unsigned int tripletIndex = rangesInGPU.tripletModuleIndices[lowerModuleIdx] + idx; + float pt = __H2F(tripletsInGPU.circleRadius[tripletIndex]) * 2 * SDL::k2Rinv1GeVf; + ; + std::vector hits = getHitsFromT3(event, tripletIndex); + SDLMath::Hit hitA(hitsInGPU.xs[hits[0]], hitsInGPU.ys[hits[0]], hitsInGPU.zs[hits[0]]); + SDLMath::Hit hitB(hitsInGPU.xs[hits[hits.size() - 1]], + hitsInGPU.ys[hits[hits.size() - 1]], + hitsInGPU.zs[hits[hits.size() - 1]]); + float eta = hitB.eta(); + float phi = hitA.phi(); + + std::vector hit_idx; + std::vector hit_type; + std::tie(hit_idx, hit_type) = getHitIdxsAndHitTypesFromT3(event, tripletIndex); + std::vector module_idx = getModuleIdxsFromT3(event, tripletIndex); + + int layer_binary = 0; + int moduleType_binary = 0; + for (size_t i = 0; i < module_idx.size(); i += 2) { + layer_binary |= (1 << (modulesInGPU.layers[module_idx[i]] + 6 * (modulesInGPU.subdets[module_idx[i]] == 4))); + moduleType_binary |= (modulesInGPU.moduleType[module_idx[i]] << i); + } + + std::vector matchedfracs; + std::vector simidx = + matchedSimTrkIdxs(hit_idx, hit_type, ana.verbose, ana.fmatch_threshold, true, matchedfracs); + + ana.tx->pushbackToBranch("t3_isFake", static_cast(simidx.size() == 0)); + ana.tx->pushbackToBranch("t3_pt", pt); + ana.tx->pushbackToBranch("t3_eta", eta); + ana.tx->pushbackToBranch("t3_phi", phi); + ana.tx->pushbackToBranch("t3_circleRadius", __H2F(tripletsInGPU.circleRadius[tripletIndex])); + ana.tx->pushbackToBranch("t3_layer_binary", layer_binary); + ana.tx->pushbackToBranch("t3_moduleType_binary", moduleType_binary); + ana.tx->pushbackToBranch("t3_partOfPT3", __H2F(tripletsInGPU.partOfPT3[tripletIndex])); + ana.tx->pushbackToBranch("t3_partOfPT5", __H2F(tripletsInGPU.partOfPT5[tripletIndex])); + ana.tx->pushbackToBranch("t3_partOfT5", __H2F(tripletsInGPU.partOfT5[tripletIndex])); + + t3_matched_simIdx.push_back(simidx); + } + } + + std::set accepted_simidxs; + for (unsigned int i = 0; i < t3_matched_simIdx.size(); ++i) { + for (unsigned int isim = 0; isim < t3_matched_simIdx[i].size(); ++isim) { + accepted_simidxs.insert(t3_matched_simIdx[i][isim]); + } + } + int n_accepted_simtrk = accepted_simidxs.size(); + std::vector sim_t3_matched(n_accepted_simtrk); + for (size_t i = 0; i < t3_matched_simIdx.size(); ++i) { + for (auto& simtrk : t3_matched_simIdx[i]) { + if (simtrk < n_accepted_simtrk) { + sim_t3_matched.at(simtrk) += 1; + } + } + } + + vector t3_isDuplicate(t3_matched_simIdx.size()); + for (unsigned int i = 0; i < t3_matched_simIdx.size(); i++) { + bool isDuplicate = false; + for (unsigned int isim = 0; isim < t3_matched_simIdx[i].size(); isim++) { + int simidx = t3_matched_simIdx[i][isim]; + if (simidx < n_accepted_simtrk) { + if (sim_t3_matched[simidx] > 1) { + isDuplicate = true; + } + } + } + t3_isDuplicate[i] = isDuplicate; + } + ana.tx->setBranch>("sim_T3_matched", sim_t3_matched); + ana.tx->setBranch>>("t3_matched_simIdx", t3_matched_simIdx); + ana.tx->setBranch>("t3_isDuplicate", t3_isDuplicate); +} + +//________________________________________________________________________________________________________________________________ +void setLowLevelBranches(SDL::Event* event) { // Get relevant information SDL::segmentsBuffer& segmentsInGPU = (*event->getSegments()); SDL::miniDoubletsBuffer& miniDoubletsInGPU = (*event->getMiniDoublets()); @@ -608,7 +697,7 @@ void setGnnNtupleBranches(SDL::Event* event) { std::vector hitidxs; std::vector hittypes; std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromTC(event, idx); - std::vector simidxs = matchedSimTrkIdxs(hitidxs, hittypes); + std::vector simidxs = matchedSimTrkIdxs(hitidxs, hittypes, ana.verbose, ana.fmatch_threshold); if (simidxs.size() == 0) continue; @@ -625,6 +714,7 @@ void setGnnNtupleBranches(SDL::Event* event) { // std::cout << " nTotalMD: " << nTotalMD << std::endl; // std::cout << " nTotalLS: " << nTotalLS << std::endl; + std::vector> LS_matched_simIdx; // Loop over modules (lower ones where the MDs are saved) for (unsigned int idx = 0; idx < *(modulesInGPU.nLowerModules); ++idx) { // // Loop over minidoublets @@ -633,7 +723,7 @@ void setGnnNtupleBranches(SDL::Event* event) { // // Get the actual index to the mini-doublet using rangesInGPU // unsigned int mdIdx = rangesInGPU.miniDoubletModuleIndices[idx] + jdx; - // setGnnNtupleMiniDoublet(event, mdIdx); + // setLowLevelMiniDoublet(event, mdIdx); // } // Loop over segments @@ -647,13 +737,13 @@ void setGnnNtupleBranches(SDL::Event* event) { if (mds_used_in_sg.find(MDs[0]) == mds_used_in_sg.end()) { mds_used_in_sg.insert(MDs[0]); md_index_map[MDs[0]] = mds_used_in_sg.size() - 1; - setGnnNtupleMiniDoublet(event, MDs[0]); + setLowLevelMiniDoublet(event, MDs[0]); } if (mds_used_in_sg.find(MDs[1]) == mds_used_in_sg.end()) { mds_used_in_sg.insert(MDs[1]); md_index_map[MDs[1]] = mds_used_in_sg.size() - 1; - setGnnNtupleMiniDoublet(event, MDs[1]); + setLowLevelMiniDoublet(event, MDs[1]); } ana.tx->pushbackToBranch("LS_MD_idx0", md_index_map[MDs[0]]); @@ -679,24 +769,15 @@ void setGnnNtupleBranches(SDL::Event* event) { std::vector hitidxs; std::vector hittypes; std::tie(hitidxs, hittypes) = getHitIdxsAndHitTypesFromLS(event, sgIdx); - std::vector simidxs = matchedSimTrkIdxs(hitidxs, hittypes); + std::vector matchedfracs; + // For LSs, set min hit matching fraction to 0 (R&D - July 2024), and write hit matching fractions + std::vector simidxs = matchedSimTrkIdxs(hitidxs, hittypes, ana.verbose, 0.0, true, matchedfracs); + + LS_matched_simIdx.push_back(simidxs); ana.tx->pushbackToBranch("LS_isFake", simidxs.size() == 0); - ana.tx->pushbackToBranch("LS_sim_pt", simidxs.size() > 0 ? trk.sim_pt()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_eta", simidxs.size() > 0 ? trk.sim_eta()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_phi", simidxs.size() > 0 ? trk.sim_phi()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pca_dxy", simidxs.size() > 0 ? trk.sim_pca_dxy()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pca_dz", simidxs.size() > 0 ? trk.sim_pca_dz()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_q", simidxs.size() > 0 ? trk.sim_q()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_event", simidxs.size() > 0 ? trk.sim_event()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_bx", simidxs.size() > 0 ? trk.sim_bunchCrossing()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_pdgId", simidxs.size() > 0 ? trk.sim_pdgId()[simidxs[0]] : -999); - ana.tx->pushbackToBranch("LS_sim_vx", - simidxs.size() > 0 ? trk.simvtx_x()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999); - ana.tx->pushbackToBranch("LS_sim_vy", - simidxs.size() > 0 ? trk.simvtx_y()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999); - ana.tx->pushbackToBranch("LS_sim_vz", - simidxs.size() > 0 ? trk.simvtx_z()[trk.sim_parentVtxIdx()[simidxs[0]]] : -999); + ana.tx->pushbackToBranch("LS_nsimMatch", simidxs.size()); + ana.tx->pushbackToBranch>("LS_matched_frac", matchedfracs); ana.tx->pushbackToBranch("LS_isInTrueTC", lss_used_in_true_tc.find(sgIdx) != lss_used_in_true_tc.end()); sg_index_map[sgIdx] = ana.tx->getBranch>("LS_isFake").size() - 1; @@ -708,6 +789,39 @@ void setGnnNtupleBranches(SDL::Event* event) { } } + std::set accepted_simidxs; + for (unsigned int i = 0; i < LS_matched_simIdx.size(); ++i) { + for (unsigned int isim = 0; isim < LS_matched_simIdx[i].size(); ++isim) { + accepted_simidxs.insert(LS_matched_simIdx[i][isim]); + } + } + int n_accepted_simtrk = accepted_simidxs.size(); + std::vector sim_LS_matched(n_accepted_simtrk); + for (size_t i = 0; i < LS_matched_simIdx.size(); ++i) { + for (auto& simtrk : LS_matched_simIdx[i]) { + if (simtrk < n_accepted_simtrk) { + sim_LS_matched.at(simtrk) += 1; + } + } + } + + vector LS_isDuplicate(LS_matched_simIdx.size()); + for (unsigned int i = 0; i < LS_matched_simIdx.size(); i++) { + bool isDuplicate = false; + for (unsigned int isim = 0; isim < LS_matched_simIdx[i].size(); isim++) { + int simidx = LS_matched_simIdx[i][isim]; + if (simidx < n_accepted_simtrk) { + if (sim_LS_matched[simidx] > 1) { + isDuplicate = true; + } + } + } + LS_isDuplicate[i] = isDuplicate; + } + ana.tx->setBranch>("sim_LS_matched", sim_LS_matched); + ana.tx->setBranch>>("LS_matched_simIdx", LS_matched_simIdx); + ana.tx->setBranch>("LS_isDuplicate", LS_isDuplicate); + for (unsigned int idx = 0; idx < nTrackCandidates; idx++) { std::vector LSs = getLSsFromTC(event, idx); std::vector lsIdx; @@ -721,7 +835,7 @@ void setGnnNtupleBranches(SDL::Event* event) { } //________________________________________________________________________________________________________________________________ -void setGnnNtupleMiniDoublet(SDL::Event* event, unsigned int MD) { +void setLowLevelMiniDoublet(SDL::Event* event, unsigned int MD) { // Get relevant information SDL::miniDoubletsBuffer& miniDoubletsInGPU = (*event->getMiniDoublets()); SDL::hitsBuffer& hitsInGPU = (*event->getHits()); @@ -743,7 +857,7 @@ void setGnnNtupleMiniDoublet(SDL::Event* event, unsigned int MD) { // Do sim matching std::vector hit_idx = {hitsInGPU.idxs[hit0], hitsInGPU.idxs[hit1]}; std::vector hit_type = {4, 4}; - std::vector simidxs = matchedSimTrkIdxs(hit_idx, hit_type); + std::vector simidxs = matchedSimTrkIdxs(hit_idx, hit_type, ana.verbose, ana.fmatch_threshold); bool isFake = simidxs.size() == 0; int tp_type = getDenomSimTrkType(simidxs); @@ -816,7 +930,7 @@ std::tuple> parseTrackCandidate(SDL:: } // Perform matching - std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type); + std::vector simidx = matchedSimTrkIdxs(hit_idx, hit_type, ana.verbose, ana.fmatch_threshold); int isFake = simidx.size() == 0; return {type, pt, eta, phi, isFake, simidx}; diff --git a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.h b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.h index aed30d5e9ebaa..2e8ebeb1e69f5 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.h +++ b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.h @@ -15,7 +15,7 @@ void createOutputBranches(); void createRequiredOutputBranches(); void createOptionalOutputBranches(); -void createGnnNtupleBranches(); +void createLowLevelBranches(); void fillOutputBranches(SDL::Event* event); void setOutputBranches(SDL::Event* event); @@ -24,8 +24,9 @@ void setOccupancyBranches(SDL::Event* event); void setPixelQuintupletOutputBranches(SDL::Event* event); void setQuintupletOutputBranches(SDL::Event* event); void setPixelTripletOutputBranches(SDL::Event* event); -void setGnnNtupleBranches(SDL::Event* event); -void setGnnNtupleMiniDoublet(SDL::Event* event, unsigned int MD); +void setTripletOutputBranches(SDL::Event* event); +void setLowLevelBranches(SDL::Event* event); +void setLowLevelMiniDoublet(SDL::Event* event, unsigned int MD); std::tuple> parseTrackCandidate(SDL::Event* event, unsigned int); std::tuple, vector> parsepT5(SDL::Event* event,