From 05289f8a678943eeab397a0bcfd57c040beccffa Mon Sep 17 00:00:00 2001 From: Chris Hoang Date: Mon, 8 Jul 2024 02:59:35 -0700 Subject: [PATCH 1/2] Address comments --- .../LSTCore/standalone/code/core/trkCore.cc | 357 +++++++++++------- .../LSTCore/standalone/code/core/trkCore.h | 7 +- .../standalone/code/core/write_sdl_ntuple.cc | 33 +- 3 files changed, 246 insertions(+), 151 deletions(-) diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc index 1a8593d12d113..a39a09a34adb8 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc @@ -274,167 +274,256 @@ float runTrackCandidate(SDL::Event *event, bool no_pls_dupclean, bool // ---------------------------------- =========================================== ---------------------------------------------- //___________________________________________________________________________________________________________________________________________________________________________________________ -std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose) { - std::vector hitidxs_(std::begin(hitidxs), std::end(hitidxs)); - std::vector hittypes_(std::begin(hittypes), std::end(hittypes)); - return matchedSimTrkIdxs(hitidxs_, hittypes_, verbose); +std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) { + std::vector hitidxs_(std::begin(hitidxs), std::end(hitidxs)); + std::vector hittypes_(std::begin(hittypes), std::end(hittypes)); + return matchedSimTrkIdxs(hitidxs_, hittypes_, verbose, matchfrac); } + //___________________________________________________________________________________________________________________________________________________________________________________________ -std::vector matchedSimTrkIdxs(std::vector hitidxs, - std::vector hittypes, - bool verbose) { - 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; - std::cout << "hittypes.size(): " << hittypes.size() << std::endl; - } +std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) +{ + std::vector matched_idxs; + std::vector matched_idx_fracs; + std::tie(matched_idxs, matched_idx_fracs) = matchedSimTrkIdxsAndFracs(hitidxs, hittypes, verbose, matchfrac); + return matched_idxs; +} - std::vector> to_check_duplicate; - for (size_t i = 0; i < hitidxs.size(); ++i) { - auto hitidx = hitidxs[i]; - auto hittype = hittypes[i]; - auto item = std::make_pair(hitidx, hittype); - if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) { - to_check_duplicate.push_back(item); +//___________________________________________________________________________________________________________________________________________________________________________________________ +std::tuple, std::vector> matchedSimTrkIdxsAndFracs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) +{ + 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; + std::cout << "hittypes.size(): " << hittypes.size() << std::endl; } - } - - int nhits_input = to_check_duplicate.size(); - std::vector> simtrk_idxs; - std::vector unique_idxs; // to aggregate which ones to count and test + std::vector> to_check_duplicate; + for (int i = 0; i < hitidxs.size(); i++) + { + auto item = std::make_pair(hitidxs[i], hittypes[i]); + if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) + { + to_check_duplicate.push_back(item); + } + } - if (verbose) { - std::cout << " '------------------------': " - << "------------------------" << std::endl; - } + int nhits_input = to_check_duplicate.size(); - for (size_t ihit = 0; ihit < to_check_duplicate.size(); ++ihit) { - auto ihitdata = to_check_duplicate[ihit]; - auto &&[hitidx, hittype] = ihitdata; + std::vector> simtrk_idxs; + std::vector unique_idxs; // to aggregate which ones to count and test - if (verbose) { - std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl; + if (verbose) + { + std::cout << " '------------------------': " + << "------------------------" << std::endl; } - std::vector simtrk_idxs_per_hit; + for (int i = 0; i < nhits_input; i++) + { + auto hitidx = to_check_duplicate[i].first; + auto hittype = to_check_duplicate[i].second; - const std::vector> *simHitIdxs = hittype == 4 ? &trk.ph2_simHitIdx() : &trk.pix_simHitIdx(); + if (verbose) + { + std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl; + } - if (verbose) { - std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; - std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; - } + std::vector simtrk_idxs_per_hit; - if (static_cast((*simHitIdxs).size()) <= hitidx) { - std::cout << "ERROR" << std::endl; - std::cout << " hittype: " << hittype << std::endl; - std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; - std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; - std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; - std::cout << hitidx << " " << hittype << std::endl; - } + const std::vector> *simHitIdxs = hittype == 4 ? &trk.ph2_simHitIdx() : &trk.pix_simHitIdx(); - for (auto &simhit_idx : (*simHitIdxs).at(hitidx)) { - if (static_cast(trk.simhit_simTrkIdx().size()) <= simhit_idx) { - std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; - std::cout << hitidx << " " << hittype << std::endl; - std::cout << trk.simhit_simTrkIdx().size() << " " << simhit_idx << std::endl; - } - int simtrk_idx = trk.simhit_simTrkIdx().at(simhit_idx); - if (verbose) { - std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx - << std::endl; - } - simtrk_idxs_per_hit.push_back(simtrk_idx); - if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end()) - unique_idxs.push_back(simtrk_idx); - } + if (verbose) + { + std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; + std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; + } - if (simtrk_idxs_per_hit.size() == 0) { - if (verbose) { - std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl; - } - simtrk_idxs_per_hit.push_back(-1); - if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end()) - unique_idxs.push_back(-1); - } + if (static_cast((*simHitIdxs).size()) <= hitidx) + { + std::cout << "ERROR" << std::endl; + std::cout << " hittype: " << hittype << std::endl; + std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; + std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; + std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; + std::cout << hitidx << " " << hittype << std::endl; + } - simtrk_idxs.push_back(simtrk_idxs_per_hit); - } + for (auto &simhit_idx : (*simHitIdxs).at(hitidx)) + { + if (static_cast(trk.simhit_simTrkIdx().size()) <= simhit_idx) + { + std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; + std::cout << hitidx << " " << hittype << std::endl; + std::cout << trk.simhit_simTrkIdx().size() << " " << simhit_idx << std::endl; + } + int simtrk_idx = trk.simhit_simTrkIdx().at(simhit_idx); + if (verbose) + { + std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx << std::endl; + } + simtrk_idxs_per_hit.push_back(simtrk_idx); + if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end()) unique_idxs.push_back(simtrk_idx); + } - if (verbose) { - std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl; - for (auto &unique_idx : unique_idxs) { - std::cout << " unique_idx: " << unique_idx << std::endl; + if (simtrk_idxs_per_hit.size() == 0) + { + if (verbose) + { + std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl; + } + simtrk_idxs_per_hit.push_back(-1); + if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end()) unique_idxs.push_back(-1); + } + + simtrk_idxs.push_back(simtrk_idxs_per_hit); } - } - // print - if (verbose) { - std::cout << "va print" << std::endl; - for (auto &vec : simtrk_idxs) { - for (auto &idx : vec) { - std::cout << idx << " "; - } - std::cout << std::endl; + if (verbose) + { + std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl; + for (auto &unique_idx : unique_idxs) + { + std::cout << " unique_idx: " << unique_idx << std::endl; + } + } + + // print + if (verbose) + { + std::cout << "va print" << std::endl; + for (auto &vec : simtrk_idxs) + { + for (auto &idx : vec) + { + std::cout << idx << " "; + } + std::cout << std::endl; + } + std::cout << "va print end" << std::endl; } - std::cout << "va print end" << std::endl; - } - // Compute all permutations - std::function> &, vector, size_t, vector> &)> perm = - [&](vector> &result, vector intermediate, size_t n, vector> &va) { + // Compute all permutations + std::function> &, vector, size_t, vector> &)> perm = [&](vector> &result, vector intermediate, size_t n, vector> &va) { // std::cout << " 'called': " << "called" << std::endl; - if (va.size() > n) { - for (auto x : va[n]) { - // std::cout << " n: " << n << std::endl; - // std::cout << " intermediate.size(): " << intermediate.size() << std::endl; - std::vector copy_intermediate(intermediate); - copy_intermediate.push_back(x); - perm(result, copy_intermediate, n + 1, va); - } - } else { - result.push_back(intermediate); + if (va.size() > n) + { + for (auto x : va[n]) + { + // std::cout << " n: " << n << std::endl; + // std::cout << " intermediate.size(): " << intermediate.size() << std::endl; + std::vector copy_intermediate(intermediate); + copy_intermediate.push_back(x); + perm(result, copy_intermediate, n + 1, va); + } + } + else + { + result.push_back(intermediate); + } + }; + + vector> allperms; + perm(allperms, vector(), 0, simtrk_idxs); + + if (verbose) + { + std::cout << " allperms.size(): " << allperms.size() << std::endl; + for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) + { + std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl; + for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) + { + std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl; + } + } + } + int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits + std::vector matched_sim_trk_idxs; + std::vector matched_sim_trk_idxs_frac; + for (auto &trkidx_perm : allperms) + { + std::vector counts; + for (auto &unique_idx : unique_idxs) + { + int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx); + counts.push_back(cnt); + } + auto result = std::max_element(counts.begin(), counts.end()); + int rawidx = std::distance(counts.begin(), result); + int trkidx = unique_idxs[rawidx]; + if (trkidx < 0) continue; + // std::cout << " counts[rawidx]: " << counts[rawidx] << std::endl; + // std::cout << " nhits_input: " << nhits_input << std::endl; + if (verbose) + { + float fr = ((float)counts[rawidx]) / ((float)nhits_input); + std::cout << " fr: " << fr << std::endl; } - }; + // std::cout << " matchfrac: " << matchfrac << std::endl; + if (counts[rawidx] > (((float)nhits_input) * matchfrac)) + { + matched_sim_trk_idxs.push_back(trkidx); + matched_sim_trk_idxs_frac.push_back( ((float)counts[rawidx]) / ((float)nhits_input) ); + } + maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end())); + } + map pairs; + unsigned size = matched_sim_trk_idxs.size(); + for (unsigned i = 0; i < size; ++i) + { + int idx = matched_sim_trk_idxs[i]; + float frac = matched_sim_trk_idxs_frac[i]; + if (pairs.find(idx) != pairs.end()) + { + if (pairs[idx] < frac) + pairs[idx] = frac; + } + else + { + pairs[idx] = frac; + } + // if (std::find(result.begin(), result.end(), idx) == result.end()) + // { + // result.push_back(idx); + // result_frac.push_back(frac); + // } + } - vector> allperms; - perm(allperms, vector(), 0, simtrk_idxs); + vector result; + vector result_frac; - if (verbose) { - std::cout << " allperms.size(): " << allperms.size() << std::endl; - for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) { - std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl; - for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) { - std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl; - } + // Loop over the map using range-based for loop + for (const auto& pair : pairs) + { + result.push_back(pair.first); + result_frac.push_back(pair.second); } - } - int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits - std::vector matched_sim_trk_idxs; - for (auto &trkidx_perm : allperms) { - std::vector counts; - for (auto &unique_idx : unique_idxs) { - int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx); - counts.push_back(cnt); - } - auto result = std::max_element(counts.begin(), counts.end()); - int rawidx = std::distance(counts.begin(), result); - int trkidx = unique_idxs[rawidx]; - if (trkidx < 0) - continue; - if (counts[rawidx] > (((float)nhits_input) * 0.75)) - matched_sim_trk_idxs.push_back(trkidx); - 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()); - return matched_sim_trk_idxs; + + // Sort indices based on 'values' + // auto indices = sort_indices(result_frac); + + // // Reorder 'vec1' and 'vec2' based on the sorted indices + // std::vector sorted_result(result.size()); + // std::vector sorted_result_frac(result_frac.size()); + // for (size_t i = 0; i < indices.size(); ++i) + // { + // sorted_result[i] = result[indices[i]]; + // sorted_result_frac[i] = result_frac[indices[i]]; + // } + + // for (size_t i = 0; i < sorted_result.size(); ++i) + // { + // int sorted_idx = sorted_result.at(i); + // float sorted_frac = sorted_result_frac.at(i); + // std::cout << " i: " << i << " sorted_idx: " << sorted_idx << " sorted_frac: " << sorted_frac << " hitidxs.size(): " << hitidxs.size() << std::endl; + // } + + return std::make_tuple(result, result_frac); + // matched_sim_trk_idxs.assign(s.begin(), s.end()); + // 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..9ee01c575bf0c 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.h +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.h @@ -28,10 +28,9 @@ float runpT3(SDL::Event *event); // --------------------- ======================== --------------------- -std::vector matchedSimTrkIdxs(std::vector hitidxs, - std::vector hittypes, - bool verbose = false); -std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose = false); +std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose=false, float matchfrac=0.75); +std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose=false, float matchfrac=0.75); +std::tuple, std::vector> matchedSimTrkIdxsAndFracs(std::vector hitidxs, std::vector hittypes, bool verbose=false, float matchfrac=0.75); 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 d8e4ce6c38ef6..be71e5219446e 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc @@ -133,7 +133,8 @@ void createOptionalOutputBranches() { ana.tx->createBranch>("t5_nonAnchorChiSquared"); // T3 branches - ana.tx->createBranch>("T3_residual"); + ana.tx->createBranch>("t3_residual"); + ana.tx->createBranch>("t3_isFake"); // Occupancy branches ana.tx->createBranch>("module_layers"); @@ -520,19 +521,25 @@ void setQuintupletOutputBranches(SDL::Event* event) { //________________________________________________________________________________________________________________________________ void setTripletOutputBranches(SDL::Event* event) { - SDL::tripletsBuffer& tripletsInGPU = (*event->getTriplets()); - SDL::objectRangesBuffer& rangesInGPU = (*event->getRanges()); - SDL::modulesBuffer& modulesInGPU = (*event->getModules()); - for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < *(modulesInGPU.nLowerModules); ++lowerModuleIdx) - { - unsigned int nTriplets = tripletsInGPU.nTriplets[lowerModuleIdx]; - for (unsigned int idx = 0; idx < nTriplets; idx++) - { - unsigned int tripletIndex = rangesInGPU.tripletModuleIndices[lowerModuleIdx] + idx; - const float residual = tripletsInGPU.residual[tripletIndex]; - ana.tx->pushbackToBranch("T3_residual", residual); - } + SDL::tripletsBuffer& tripletsInGPU = (*event->getTriplets()); + SDL::objectRangesBuffer& rangesInGPU = (*event->getRanges()); + SDL::modulesBuffer& modulesInGPU = (*event->getModules()); + for (unsigned int lowerModuleIdx = 0; lowerModuleIdx < *(modulesInGPU.nLowerModules); ++lowerModuleIdx) { + unsigned int nTriplets = tripletsInGPU.nTriplets[lowerModuleIdx]; + for (unsigned int idx = 0; idx < nTriplets; idx++) { + unsigned int tripletIndex = rangesInGPU.tripletModuleIndices[lowerModuleIdx] + idx; + + const float residual = tripletsInGPU.residual[tripletIndex]; + ana.tx->pushbackToBranch("t3_residual", residual); + + std::tuple, std::vector> hitIdxsAndTypes = getHitIdxsAndHitTypesFromT3(event, tripletIndex); + std::vector hit_idxs = std::get<0>(hitIdxsAndTypes); + std::vector hit_types = std::get<1>(hitIdxsAndTypes); + + std::tuple, std::vector> simIdxsAndFracs = matchedSimTrkIdxsAndFracs(hit_idxs, hit_types); + ana.tx->pushbackToBranch("t3_isFake", static_cast(std::get<0>(simIdxsAndFracs).size() == 0)); } + } } //________________________________________________________________________________________________________________________________ From 7ff7999787418354bd2c1dfca25a3dfae3606fcd Mon Sep 17 00:00:00 2001 From: Chris Hoang Date: Mon, 8 Jul 2024 12:47:44 -0700 Subject: [PATCH 2/2] Fix styling --- .../LSTCore/standalone/code/core/trkCore.cc | 393 ++++++++---------- 1 file changed, 179 insertions(+), 214 deletions(-) diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc index a39a09a34adb8..e215e330d5b11 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc @@ -275,255 +275,220 @@ float runTrackCandidate(SDL::Event *event, bool no_pls_dupclean, bool //___________________________________________________________________________________________________________________________________________________________________________________________ std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) { - std::vector hitidxs_(std::begin(hitidxs), std::end(hitidxs)); - std::vector hittypes_(std::begin(hittypes), std::end(hittypes)); - return matchedSimTrkIdxs(hitidxs_, hittypes_, verbose, matchfrac); + std::vector hitidxs_(std::begin(hitidxs), std::end(hitidxs)); + std::vector hittypes_(std::begin(hittypes), std::end(hittypes)); + return matchedSimTrkIdxs(hitidxs_, hittypes_, verbose, matchfrac); } //___________________________________________________________________________________________________________________________________________________________________________________________ -std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) -{ - std::vector matched_idxs; - std::vector matched_idx_fracs; - std::tie(matched_idxs, matched_idx_fracs) = matchedSimTrkIdxsAndFracs(hitidxs, hittypes, verbose, matchfrac); - return matched_idxs; +std::vector matchedSimTrkIdxs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) { + std::vector matched_idxs; + std::vector matched_idx_fracs; + std::tie(matched_idxs, matched_idx_fracs) = matchedSimTrkIdxsAndFracs(hitidxs, hittypes, verbose, matchfrac); + return matched_idxs; } //___________________________________________________________________________________________________________________________________________________________________________________________ -std::tuple, std::vector> matchedSimTrkIdxsAndFracs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) -{ - 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; - std::cout << "hittypes.size(): " << hittypes.size() << std::endl; - } +std::tuple, std::vector> matchedSimTrkIdxsAndFracs(std::vector hitidxs, std::vector hittypes, bool verbose, float matchfrac) { + 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; + std::cout << "hittypes.size(): " << hittypes.size() << std::endl; + } - std::vector> to_check_duplicate; - for (int i = 0; i < hitidxs.size(); i++) - { - auto item = std::make_pair(hitidxs[i], hittypes[i]); - if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) - { - to_check_duplicate.push_back(item); - } + std::vector> to_check_duplicate; + for (int i = 0; i < hitidxs.size(); i++) { + auto item = std::make_pair(hitidxs[i], hittypes[i]); + if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) { + to_check_duplicate.push_back(item); } + } - int nhits_input = to_check_duplicate.size(); - - std::vector> simtrk_idxs; - std::vector unique_idxs; // to aggregate which ones to count and test - - if (verbose) - { - std::cout << " '------------------------': " - << "------------------------" << std::endl; - } + int nhits_input = to_check_duplicate.size(); - for (int i = 0; i < nhits_input; i++) - { - auto hitidx = to_check_duplicate[i].first; - auto hittype = to_check_duplicate[i].second; + std::vector> simtrk_idxs; + std::vector unique_idxs; // to aggregate which ones to count and test - if (verbose) - { - std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl; - } + if (verbose) { + std::cout << " '------------------------': " + << "------------------------" << std::endl; + } - std::vector simtrk_idxs_per_hit; + for (int i = 0; i < nhits_input; i++) { + auto hitidx = to_check_duplicate[i].first; + auto hittype = to_check_duplicate[i].second; - const std::vector> *simHitIdxs = hittype == 4 ? &trk.ph2_simHitIdx() : &trk.pix_simHitIdx(); + if (verbose) { + std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl; + } - if (verbose) - { - std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; - std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; - } + std::vector simtrk_idxs_per_hit; - if (static_cast((*simHitIdxs).size()) <= hitidx) - { - std::cout << "ERROR" << std::endl; - std::cout << " hittype: " << hittype << std::endl; - std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; - std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; - std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; - std::cout << hitidx << " " << hittype << std::endl; - } + const std::vector> *simHitIdxs = hittype == 4 ? &trk.ph2_simHitIdx() : &trk.pix_simHitIdx(); - for (auto &simhit_idx : (*simHitIdxs).at(hitidx)) - { - if (static_cast(trk.simhit_simTrkIdx().size()) <= simhit_idx) - { - std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; - std::cout << hitidx << " " << hittype << std::endl; - std::cout << trk.simhit_simTrkIdx().size() << " " << simhit_idx << std::endl; - } - int simtrk_idx = trk.simhit_simTrkIdx().at(simhit_idx); - if (verbose) - { - std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx << std::endl; - } - simtrk_idxs_per_hit.push_back(simtrk_idx); - if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end()) unique_idxs.push_back(simtrk_idx); - } + if (verbose) { + std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; + std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; + } - if (simtrk_idxs_per_hit.size() == 0) - { - if (verbose) - { - std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl; - } - simtrk_idxs_per_hit.push_back(-1); - if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end()) unique_idxs.push_back(-1); - } + if (static_cast((*simHitIdxs).size()) <= hitidx) { + std::cout << "ERROR" << std::endl; + std::cout << " hittype: " << hittype << std::endl; + std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl; + std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl; + std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; + std::cout << hitidx << " " << hittype << std::endl; + } - simtrk_idxs.push_back(simtrk_idxs_per_hit); + for (auto &simhit_idx : (*simHitIdxs).at(hitidx)) { + if (static_cast(trk.simhit_simTrkIdx().size()) <= simhit_idx) { + std::cout << (*simHitIdxs).size() << " " << hittype << std::endl; + std::cout << hitidx << " " << hittype << std::endl; + std::cout << trk.simhit_simTrkIdx().size() << " " << simhit_idx << std::endl; + } + int simtrk_idx = trk.simhit_simTrkIdx().at(simhit_idx); + if (verbose) { + std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx << std::endl; + } + simtrk_idxs_per_hit.push_back(simtrk_idx); + if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end()) unique_idxs.push_back(simtrk_idx); } - if (verbose) - { - std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl; - for (auto &unique_idx : unique_idxs) - { - std::cout << " unique_idx: " << unique_idx << std::endl; - } + if (simtrk_idxs_per_hit.size() == 0) { + if (verbose) { + std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl; + } + simtrk_idxs_per_hit.push_back(-1); + if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end()) unique_idxs.push_back(-1); } + simtrk_idxs.push_back(simtrk_idxs_per_hit); + } - // print - if (verbose) - { - std::cout << "va print" << std::endl; - for (auto &vec : simtrk_idxs) - { - for (auto &idx : vec) - { - std::cout << idx << " "; - } - std::cout << std::endl; - } - std::cout << "va print end" << std::endl; + if (verbose) { + std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl; + for (auto &unique_idx : unique_idxs) { + std::cout << " unique_idx: " << unique_idx << std::endl; } + } - // Compute all permutations - std::function> &, vector, size_t, vector> &)> perm = [&](vector> &result, vector intermediate, size_t n, vector> &va) { - // std::cout << " 'called': " << "called" << std::endl; - if (va.size() > n) - { - for (auto x : va[n]) - { - // std::cout << " n: " << n << std::endl; - // std::cout << " intermediate.size(): " << intermediate.size() << std::endl; - std::vector copy_intermediate(intermediate); - copy_intermediate.push_back(x); - perm(result, copy_intermediate, n + 1, va); - } - } - else - { - result.push_back(intermediate); - } - }; - - vector> allperms; - perm(allperms, vector(), 0, simtrk_idxs); - - if (verbose) - { - std::cout << " allperms.size(): " << allperms.size() << std::endl; - for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) - { - std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl; - for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) - { - std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl; - } - } + // print + if (verbose) { + std::cout << "va print" << std::endl; + for (auto &vec : simtrk_idxs) { + for (auto &idx : vec) { + std::cout << idx << " "; + } + std::cout << std::endl; } - int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits - std::vector matched_sim_trk_idxs; - std::vector matched_sim_trk_idxs_frac; - for (auto &trkidx_perm : allperms) - { - std::vector counts; - for (auto &unique_idx : unique_idxs) - { - int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx); - counts.push_back(cnt); - } - auto result = std::max_element(counts.begin(), counts.end()); - int rawidx = std::distance(counts.begin(), result); - int trkidx = unique_idxs[rawidx]; - if (trkidx < 0) continue; - // std::cout << " counts[rawidx]: " << counts[rawidx] << std::endl; - // std::cout << " nhits_input: " << nhits_input << std::endl; - if (verbose) - { - float fr = ((float)counts[rawidx]) / ((float)nhits_input); - std::cout << " fr: " << fr << std::endl; - } - // std::cout << " matchfrac: " << matchfrac << std::endl; - if (counts[rawidx] > (((float)nhits_input) * matchfrac)) - { - matched_sim_trk_idxs.push_back(trkidx); - matched_sim_trk_idxs_frac.push_back( ((float)counts[rawidx]) / ((float)nhits_input) ); - } - maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end())); + std::cout << "va print end" << std::endl; + } + + // Compute all permutations + std::function> &, vector, size_t, vector> &)> perm = [&](vector> &result, vector intermediate, size_t n, vector> &va) { + // std::cout << " 'called': " << "called" << std::endl; + if (va.size() > n) { + for (auto x : va[n]) { + // std::cout << " n: " << n << std::endl; + // std::cout << " intermediate.size(): " << intermediate.size() << std::endl; + std::vector copy_intermediate(intermediate); + copy_intermediate.push_back(x); + perm(result, copy_intermediate, n + 1, va); + } } - map pairs; - unsigned size = matched_sim_trk_idxs.size(); - for (unsigned i = 0; i < size; ++i) - { - int idx = matched_sim_trk_idxs[i]; - float frac = matched_sim_trk_idxs_frac[i]; - if (pairs.find(idx) != pairs.end()) - { - if (pairs[idx] < frac) - pairs[idx] = frac; - } - else - { - pairs[idx] = frac; - } - // if (std::find(result.begin(), result.end(), idx) == result.end()) - // { - // result.push_back(idx); - // result_frac.push_back(frac); - // } + else { + result.push_back(intermediate); } + }; - vector result; - vector result_frac; + vector> allperms; + perm(allperms, vector(), 0, simtrk_idxs); - // Loop over the map using range-based for loop - for (const auto& pair : pairs) - { - result.push_back(pair.first); - result_frac.push_back(pair.second); + if (verbose) { + std::cout << " allperms.size(): " << allperms.size() << std::endl; + for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) { + std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl; + for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) { + std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl; + } } - - // Sort indices based on 'values' - // auto indices = sort_indices(result_frac); - - // // Reorder 'vec1' and 'vec2' based on the sorted indices - // std::vector sorted_result(result.size()); - // std::vector sorted_result_frac(result_frac.size()); - // for (size_t i = 0; i < indices.size(); ++i) + } + int maxHitMatchCount = 0; // ultimate maximum of the number of matched hits + std::vector matched_sim_trk_idxs; + std::vector matched_sim_trk_idxs_frac; + for (auto &trkidx_perm : allperms) { + std::vector counts; + for (auto &unique_idx : unique_idxs) { + int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx); + counts.push_back(cnt); + } + auto result = std::max_element(counts.begin(), counts.end()); + int rawidx = std::distance(counts.begin(), result); + int trkidx = unique_idxs[rawidx]; + if (trkidx < 0) continue; + // std::cout << " counts[rawidx]: " << counts[rawidx] << std::endl; + // std::cout << " nhits_input: " << nhits_input << std::endl; + if (verbose) { + float fr = ((float)counts[rawidx]) / ((float)nhits_input); + std::cout << " fr: " << fr << std::endl; + } + // std::cout << " matchfrac: " << matchfrac << std::endl; + if (counts[rawidx] > (((float)nhits_input) * matchfrac)) { + matched_sim_trk_idxs.push_back(trkidx); + matched_sim_trk_idxs_frac.push_back( ((float)counts[rawidx]) / ((float)nhits_input) ); + } + maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end())); + } + map pairs; + unsigned size = matched_sim_trk_idxs.size(); + for (unsigned i = 0; i < size; ++i) { + int idx = matched_sim_trk_idxs[i]; + float frac = matched_sim_trk_idxs_frac[i]; + if (pairs.find(idx) != pairs.end()) { + if (pairs[idx] < frac) + pairs[idx] = frac; + } + else { + pairs[idx] = frac; + } + // if (std::find(result.begin(), result.end(), idx) == result.end()) // { - // sorted_result[i] = result[indices[i]]; - // sorted_result_frac[i] = result_frac[indices[i]]; + // result.push_back(idx); + // result_frac.push_back(frac); // } + } - // for (size_t i = 0; i < sorted_result.size(); ++i) - // { - // int sorted_idx = sorted_result.at(i); - // float sorted_frac = sorted_result_frac.at(i); - // std::cout << " i: " << i << " sorted_idx: " << sorted_idx << " sorted_frac: " << sorted_frac << " hitidxs.size(): " << hitidxs.size() << std::endl; - // } + vector result; + vector result_frac; + + // Loop over the map using range-based for loop + for (const auto& pair : pairs) { + result.push_back(pair.first); + result_frac.push_back(pair.second); + } - return std::make_tuple(result, result_frac); - // matched_sim_trk_idxs.assign(s.begin(), s.end()); - // return matched_sim_trk_idxs; + // Sort indices based on 'values' + // auto indices = sort_indices(result_frac); + + // // Reorder 'vec1' and 'vec2' based on the sorted indices + // std::vector sorted_result(result.size()); + // std::vector sorted_result_frac(result_frac.size()); + // for (size_t i = 0; i < indices.size(); ++i) + // { + // sorted_result[i] = result[indices[i]]; + // sorted_result_frac[i] = result_frac[indices[i]]; + // } + + // for (size_t i = 0; i < sorted_result.size(); ++i) + // { + // int sorted_idx = sorted_result.at(i); + // float sorted_frac = sorted_result_frac.at(i); + // std::cout << " i: " << i << " sorted_idx: " << sorted_idx << " sorted_frac: " << sorted_frac << " hitidxs.size(): " << hitidxs.size() << std::endl; + // } + + return std::make_tuple(result, result_frac); + // matched_sim_trk_idxs.assign(s.begin(), s.end()); + // return matched_sim_trk_idxs; } //__________________________________________________________________________________________