Add low-level information in standalone LST output ntuple

RecoTracker/LSTCore/standalone/bin/sdl.cc

@@ -51,6 +51,7 @@ int main(int argc, char **argv) {
       "o,output", "Output file name", cxxopts::value<std::string>())(
       "N,nmatch", "N match for MTV-like matching", cxxopts::value<int>()->default_value("9"))(
       "p,ptCut", "Min pT cut In GeV", cxxopts::value<float>()->default_value("0.8"))(
+      "f,fmatch", "hit matching fraction for RECO-SIM matching", cxxopts::value<float>()->default_value("0.75"))(
       "n,nevents", "N events to loop over", cxxopts::value<int>()->default_value("-1"))(
       "x,event_index", "specific event index to process", cxxopts::value<int>()->default_value("-1"))(
       "g,pdg_id", "The simhit pdgId match option", cxxopts::value<int>()->default_value("0"))(
@@ -60,7 +61,8 @@ int main(int argc, char **argv) {
       "w,write_ntuple", "Write Ntuple", cxxopts::value<int>()->default_value("1"))(
       "s,streams", "Set number of streams", cxxopts::value<int>()->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<int>())(
       "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<int>();
 
+  //_______________________________________________________________________________
+  // --fmatch
+  ana.fmatch_threshold = result["fmatch"].as<float>();
+
   //_______________________________________________________________________________
   // --nevents
   ana.n_events = result["nevents"].as<int>();
@@ -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();
     }
   }
 

RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc

@@ -139,6 +139,33 @@ std::vector<unsigned int> getHitsFromT3(SDL::Event<SDL::Acc>* event, unsigned in
   return {hits_0[0], hits_0[1], hits_1[0], hits_1[1], hits_2[0], hits_2[1]};
 }
 
+//____________________________________________________________________________________________
+std::vector<unsigned int> getHitIdxsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3) {
+  SDL::hitsBuffer<alpaka::DevCpu>& hitsInGPU = *(event->getHits());
+  std::vector<unsigned int> hits = getHitsFromT3(event, T3);
+  std::vector<unsigned int> hitidxs;
+  for (auto& hit : hits)
+    hitidxs.push_back(hitsInGPU.idxs[hit]);
+  return hitidxs;
+}
+
+//____________________________________________________________________________________________
+std::vector<unsigned int> getModuleIdxsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3) {
+  std::vector<unsigned int> hits = getHitsFromT3(event, T3);
+  std::vector<unsigned int> module_idxs;
+  SDL::hitsBuffer<alpaka::DevCpu>& hitsInGPU = *(event->getHits());
+  for (auto& hitIdx : hits) {
+    module_idxs.push_back(hitsInGPU.moduleIndices[hitIdx]);
+  }
+  return module_idxs;
+}
+
+//____________________________________________________________________________________________
+std::vector<unsigned int> getHitTypesFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3) {
+  return {4, 4, 4, 4, 4, 4};
+  ;
+}
+
 //____________________________________________________________________________________________
 std::tuple<std::vector<unsigned int>, std::vector<unsigned int>> getHitIdxsAndHitTypesFromT3(
     SDL::Event<SDL::Acc>* event, unsigned T3) {

RecoTracker/LSTCore/standalone/code/core/AccessHelper.h

@@ -33,6 +33,9 @@ std::tuple<std::vector<unsigned int>, std::vector<unsigned int>> getHitIdxsAndHi
 std::vector<unsigned int> getLSsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
 std::vector<unsigned int> getMDsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
 std::vector<unsigned int> getHitsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
+std::vector<unsigned int> getHitIdxsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
+std::vector<unsigned int> getHitTypesFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
+std::vector<unsigned int> getModuleIdxsFromT3(SDL::Event<SDL::Acc>* event, unsigned int T3);
 std::tuple<std::vector<unsigned int>, std::vector<unsigned int>> getHitIdxsAndHitTypesFromT3(
     SDL::Event<SDL::Acc>* event, unsigned T3);
 

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;
 

RecoTracker/LSTCore/standalone/code/core/trkCore.cc

@@ -283,7 +283,10 @@ std::vector<int> matchedSimTrkIdxs(std::vector<int> hitidxs, std::vector<int> hi
 //___________________________________________________________________________________________________________________________________________________________________________________________
 std::vector<int> matchedSimTrkIdxs(std::vector<unsigned int> hitidxs,
                                    std::vector<unsigned int> hittypes,
-                                   bool verbose) {
+                                   bool verbose,
+                                   float matchFraction,
+                                   bool writeMatchFraction,
+                                   std::vector<float> &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<int> matchedSimTrkIdxs(std::vector<unsigned int> hitidxs,
   }
   int maxHitMatchCount = 0;  // ultimate maximum of the number of matched hits
   std::vector<int> matched_sim_trk_idxs;
+  std::vector<std::pair<int, float>> matched_sim_idx_frac;
   for (auto &trkidx_perm : allperms) {
     std::vector<int> counts;
     for (auto &unique_idx : unique_idxs) {
@@ -425,15 +429,22 @@ std::vector<int> matchedSimTrkIdxs(std::vector<unsigned int> 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<int> 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;
 }
 

RecoTracker/LSTCore/standalone/code/core/trkCore.h

@@ -28,9 +28,13 @@ float runpT3(SDL::Event<SDL::Acc> *event);
 
 // --------------------- ======================== ---------------------
 
+static std::vector<float> empty_vector_float = std::vector<float>();
 std::vector<int> matchedSimTrkIdxs(std::vector<unsigned int> hitidxs,
                                    std::vector<unsigned int> hittypes,
-                                   bool verbose = false);
+                                   bool verbose = false,
+                                   float matchFraction = 0.75,
+                                   bool writeMatchFraction = false,
+                                   std::vector<float> &matched_fractions = empty_vector_float);
 std::vector<int> matchedSimTrkIdxs(std::vector<int> hitidxs, std::vector<int> hittypes, bool verbose = false);
 int getDenomSimTrkType(int isimtrk);
 int getDenomSimTrkType(std::vector<int> simidxs);