feat: Lepton PID

CODEOWNERS

@@ -6,13 +6,16 @@
 # @RichCap                     Richard Capobianco
 # @rtysonCLAS12                Richard Tyson
 # @Gregtom3                    Gregory Matousek
+# @MarianaT27                  Mariana Tenorio Pita
+# @tbhayward                   Timothy Hayward
 #################################################################################
 
 * @c-dilks
 src/iguana/algorithms/clas12/CalorimeterLinker/* @c-dilks
 src/iguana/algorithms/clas12/EventBuilderFilter/* @c-dilks
 src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar
 src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3
+src/iguana/algorithms/clas12/LeptonIDFilter/* @MarianaT27
 src/iguana/algorithms/clas12/MatchParticleProximity/* @c-dilks
 src/iguana/algorithms/clas12/MomentumCorrection/* @RichCap @c-dilks
 src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3

examples/meson.build

@@ -1,5 +1,5 @@
 # install config files
-example_config_files_prefix = project_etcdir / 'examples'
+example_config_files_prefix = project_etc_dir / 'examples'
 install_subdir('config', install_dir: example_config_files_prefix, strip_directory: true)
 
 # example source information

meson.build

@@ -9,6 +9,7 @@ project(
     'cpp_std': 'c++17',
     'buildtype': 'release',
     'libdir': 'lib',
+    'datadir': 'share',
     'licensedir': 'share/licenses/iguana',
     'pkgconfig.relocatable': 'true',
     'force_fallback_for': ['rcdb'],
@@ -20,6 +21,14 @@ project_description = 'Implementation Guardian of Analysis Algorithms'
 add_languages('fortran', native: false, required: get_option('bind_fortran'))
 use_chameleon = get_option('bind_fortran')
 
+# check built-in build options
+if get_option('libdir') != 'lib'
+  warning('build option "libdir" = "' + get_option('libdir') + '", which is not "lib"; if you experience any issues, please report them')
+endif
+if get_option('datadir') != 'share'
+  error('build option "datadir" must be "share", but it is set to "' + get_option('datadir') + '"')
+endif
+
 # warn that macOS is no longer tested
 if host_machine.system() == 'darwin'
   warning('''host machine system is darwin:
@@ -70,6 +79,7 @@ ROOT_dep = dependency(
     'ROOT::RIO',
     'ROOT::ROOTDataFrame',
     'ROOT::ROOTVecOps',
+    'ROOT::TMVA',
     'ROOT::TreePlayer',
   ],
 )
@@ -154,7 +164,8 @@ hipo_dep_dataframes_found = hipo_dep.get_variable(
 project_inc = include_directories('src')
 project_libs = []
 project_deps = declare_dependency(dependencies: dep_list)
-project_etcdir = get_option('sysconfdir') / meson.project_name()
+project_etc_dir = get_option('sysconfdir') / meson.project_name()
+project_data_dir = get_option('datadir') / meson.project_name()
 project_test_env = environment()
 project_pythondir = 'python'
 
@@ -180,9 +191,14 @@ project_test_env.set(
 )
 
 # set preprocessor macros
-add_project_arguments('-DIGUANA_ETCDIR="' + get_option('prefix') / project_etcdir + '"', language: [ 'cpp' ])
+add_project_arguments(
+  '-D IGUANA_PREFIX="' + get_option('prefix') + '"',
+  '-D IGUANA_ETCDIR="' + project_etc_dir + '"',
+  '-D IGUANA_DATADIR="' + project_data_dir + '"',
+  language: ['cpp'],
+)
 if ROOT_dep.found()
-  add_project_arguments('-DIGUANA_ROOT_FOUND', language: [ 'cpp' ]) # currently only used for Validator plot styles
+  add_project_arguments('-D IGUANA_ROOT_FOUND', language: [ 'cpp' ]) # currently only used for Validator plot styles
 endif
 
 # start chameleon
@@ -236,7 +252,7 @@ if get_option('z_install_envfile')
       'libdir':            get_option('libdir'),
       'includedir':        get_option('includedir'),
       'bindir':            get_option('bindir'),
-      'etcdir':            project_etcdir,
+      'etcdir':            project_etc_dir,
       'pythondir':         project_pythondir,
       'dep_pkgconfigdirs': ':'.join(dep_pkgconfig_dirs),
       'dep_libdirs':       ':'.join(dep_lib_dirs),

meson/lsan.supp

@@ -1 +1,2 @@
 leak:libCling.so
+leak:TMVA::DataSetFactory::BuildDynamicDataSet

src/iguana/algorithms/Algorithm.cc

@@ -86,6 +86,19 @@ namespace iguana {
 
   ///////////////////////////////////////////////////////////////////////////////
 
+  YAML::Node Algorithm::GetOptionNode(YAMLReader::node_path_t node_path) const
+  {
+    CompleteOptionNodePath("", node_path);
+    auto node = m_yaml_config->GetNode(node_path);
+    if(!node.has_value()) {
+      m_log->Error("Algorithm::GetOptionNode failed to find node");
+      throw std::runtime_error("config file parsing issue");
+    }
+    return node.value();
+  }
+
+  ///////////////////////////////////////////////////////////////////////////////
+
   void Algorithm::SetName(std::string_view name)
   {
     Object::SetName(name);
@@ -123,6 +136,17 @@ namespace iguana {
 
   ///////////////////////////////////////////////////////////////////////////////
 
+  std::string Algorithm::GetDataFile(std::string const& name)
+  {
+    if(!m_datafile_reader) {
+      m_datafile_reader = std::make_unique<DataFileReader>(ConfigFileReader::ConvertAlgoNameToConfigDir(m_class_name), "data|" + m_name);
+      m_datafile_reader->SetLogLevel(m_log->GetLevel());
+    }
+    return m_datafile_reader->FindFile(name);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////////
+
   void Algorithm::ParseYAMLConfig()
   {
 

src/iguana/algorithms/Algorithm.h

@@ -8,9 +8,10 @@
 
 #include "AlgorithmBoilerplate.h"
 #include "iguana/bankdefs/BankDefs.h"
+#include "iguana/services/DataFileReader.h"
+#include "iguana/services/GlobalParam.h"
 #include "iguana/services/RCDBReader.h"
 #include "iguana/services/YAMLReader.h"
-#include <iguana/services/GlobalParam.h>
 
 namespace iguana {
 
@@ -144,6 +145,11 @@ namespace iguana {
       template <typename OPTION_TYPE>
       std::set<OPTION_TYPE> GetOptionSet(std::string const& key, YAMLReader::node_path_t node_path = {}) const;
 
+      /// Get the `YAML::Node` for an option
+      /// @param node_path the `YAML::Node` identifier path to search
+      /// @returns the `YAML::Node`
+      YAML::Node GetOptionNode(YAMLReader::node_path_t node_path) const;
+
       /// Set the name of this algorithm
       /// @param name the new name
       void SetName(std::string_view name);
@@ -166,6 +172,11 @@ namespace iguana {
       /// @param name the directory name
       void SetConfigDirectory(std::string const& name);
 
+      /// Get the full path to a data file, such as a machine-learning model
+      /// @param name the name of the file; if found in the user's current working directory (`./`), that will be the file that is used;
+      /// otherwise the _installed_ file (in `$IGUANA/share/`) will be used by default
+      std::string GetDataFile(std::string const& name);
+
       /// Get the index of a bank in a `hipo::banklist`; throws an exception if the bank is not found
       /// @param banks the list of banks this algorithm will use
       /// @param bank_name the name of the bank
@@ -311,6 +322,9 @@ namespace iguana {
 
       /// Data structure to hold configuration options set by `Algorithm::SetOption`
       std::unordered_map<std::string, option_t> m_option_cache;
+
+      /// Data file reader
+      std::unique_ptr<DataFileReader> m_datafile_reader;
   };
 
   //////////////////////////////////////////////////////////////////////////////

src/iguana/algorithms/clas12/LeptonIDFilter/Algorithm.cc

@@ -0,0 +1,207 @@
+#include "Algorithm.h"
+
+#include <Math/Vector4D.h>
+#include <cmath>
+
+namespace iguana::clas12 {
+
+  REGISTER_IGUANA_ALGORITHM(LeptonIDFilter, "clas12::LeptonIDFilter");
+
+  void LeptonIDFilter::Start(hipo::banklist& banks)
+  {
+    // Get configuration
+    ParseYAMLConfig();
+    o_pids                = GetOptionSet<int>("pids");
+    o_cut                 = GetOptionScalar<double>("cut");
+    o_tmva_reader_options = GetOptionScalar<std::string>("tmva_reader_options");
+    o_particle_bank       = GetOptionScalar<std::string>("particle_bank");
+    o_runnum              = ConcurrentParamFactory::Create<int>();
+    o_weightfile_electron = ConcurrentParamFactory::Create<std::string>();
+    o_weightfile_positron = ConcurrentParamFactory::Create<std::string>();
+
+    // Get Banks that we are going to use
+    b_particle    = GetBankIndex(banks, o_particle_bank);
+    b_calorimeter = GetBankIndex(banks, "REC::Calorimeter");
+    b_config      = GetBankIndex(banks, "RUN::config");
+
+    // Initialize the TMVA reader
+    readerTMVA = std::make_unique<TMVA::Reader>(LeptonIDVars::names, o_tmva_reader_options);
+
+    // find all the unique weights files in the configuration YAML
+    std::set<std::string> weightfile_list;
+    for(auto const& node : GetOptionNode({"weightfile"})) {
+      for(std::string const particle : {"electron", "positron"}) {
+        if(node[particle]) {
+          weightfile_list.insert(node[particle].as<std::string>());
+        }
+      }
+    }
+
+    // book the weights files
+    // use the name of weightfile as the method tag, for simplicity
+    m_log->Debug("Booking weight files:");
+    for(auto const& weightfile_name : weightfile_list) {
+      auto weightfile_path = GetDataFile(weightfile_name);
+      m_log->Debug(" - {}", weightfile_path);
+      readerTMVA->BookMVA(weightfile_name, weightfile_path);
+    }
+  }
+
+
+  bool LeptonIDFilter::Run(hipo::banklist& banks) const
+  {
+    return Run(
+        GetBank(banks, b_particle, o_particle_bank),
+        GetBank(banks, b_calorimeter, "REC::Calorimeter"),
+        GetBank(banks, b_config, "RUN::config"));
+  }
+
+
+  bool LeptonIDFilter::Run(hipo::bank& particleBank, hipo::bank const& calorimeterBank, hipo::bank const& configBank) const
+  {
+    // particle bank before filtering
+    ShowBank(particleBank, Logger::Header("INPUT PARTICLES"));
+
+    // prepare the event, reloading configuration parameters if the run number changed or is not yet known
+    auto key = PrepareEvent(configBank.getInt("run", 0));
+
+    // filter the particle bank
+    particleBank.getMutableRowList().filter([this, &calorimeterBank, key](auto bank, auto row) {
+      auto pid = bank.getInt("pid", row);
+      // check if this is a lepton in `o_pids`
+      if(o_pids.find(pid) != o_pids.end()) {
+        auto status = bank.getShort("status", row);
+        // status cut
+        if(std::abs(status) >= 2000 && std::abs(status) < 4000) {
+          m_log->Trace("Found lepton: pindex={}", row);
+          auto lepton_vars  = GetLeptonIDVariables(row, bank, calorimeterBank);
+          lepton_vars.pid   = pid;
+          lepton_vars.score = CalculateScore(lepton_vars, key);
+          return Filter(lepton_vars.score) ? 1 : 0;
+        }
+        else {
+          m_log->Trace("Lepton at pindex={} did not pass status cut", row);
+          return 0;
+        }
+      }
+      return 1; // not a lepton in `o_pids`, let it pass the filter
+    });
+
+    // particle bank after filtering
+    ShowBank(particleBank, Logger::Header("OUTPUT PARTICLES"));
+    return !particleBank.getRowList().empty();
+  }
+
+
+  concurrent_key_t LeptonIDFilter::PrepareEvent(int const runnum) const
+  {
+    m_log->Trace("calling PrepareEvent({})", runnum);
+    if(o_runnum->NeedsHashing()) {
+      std::hash<int> hash_ftn;
+      auto hash_key = hash_ftn(runnum);
+      if(!o_runnum->HasKey(hash_key))
+        Reload(runnum, hash_key);
+      return hash_key;
+    }
+    else {
+      if(o_runnum->IsEmpty() || o_runnum->Load(0) != runnum)
+        Reload(runnum, 0);
+      return 0;
+    }
+  }
+
+
+  void LeptonIDFilter::Reload(int const runnum, concurrent_key_t key) const
+  {
+    std::lock_guard<std::mutex> const lock(m_mutex); // NOTE: be sure to lock successive `ConcurrentParam::Save` calls !!!
+    m_log->Trace("-> calling Reload({}, {})", runnum, key);
+    o_runnum->Save(runnum, key);
+    o_weightfile_electron->Save(GetOptionScalar<std::string>("weightfile:electron", {"weightfile", GetConfig()->InRange("runs", runnum), "electron"}), key);
+    o_weightfile_positron->Save(GetOptionScalar<std::string>("weightfile:positron", {"weightfile", GetConfig()->InRange("runs", runnum), "positron"}), key);
+  }
+
+
+  LeptonIDVars LeptonIDFilter::GetLeptonIDVariables(int const plepton, hipo::bank const& particle_bank, hipo::bank const& calorimeter_bank) const
+  {
+
+    double px = particle_bank.getFloat("px", plepton);
+    double py = particle_bank.getFloat("py", plepton);
+    double pz = particle_bank.getFloat("pz", plepton);
+    double E  = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2) + std::pow(0.000511, 2));
+    ROOT::Math::PxPyPzMVector vec_lepton(px, py, pz, E);
+
+    LeptonIDVars lepton;
+
+    lepton.P     = vec_lepton.P();
+    lepton.Theta = vec_lepton.Theta();
+    lepton.Phi   = vec_lepton.Phi();
+
+    m_log->Debug("Variables obtained from particle bank");
+
+
+    lepton.m2pcal  = -1;
+    lepton.m2ecin  = -1;
+    lepton.m2ecout = -1;
+
+    for(int row = 0; row < calorimeter_bank.getRows(); row++) {
+      auto pindex = calorimeter_bank.getShort("pindex", row);
+      auto layer  = calorimeter_bank.getByte("layer", row);
+      auto energy = calorimeter_bank.getFloat("energy", row);
+      auto m2u    = calorimeter_bank.getFloat("m2u", row);
+      auto m2v    = calorimeter_bank.getFloat("m2v", row);
+      auto m2w    = calorimeter_bank.getFloat("m2w", row);
+
+      if(pindex == plepton && layer == 1) {
+        lepton.SFpcal = energy / vec_lepton.P();
+        lepton.m2pcal = (m2u + m2v + m2w) / 3;
+      }
+
+      if(pindex == plepton && layer == 4) {
+        lepton.SFecin = energy / vec_lepton.P();
+        lepton.m2ecin = (m2u + m2v + m2w) / 3;
+      }
+      if(pindex == plepton && layer == 7) {
+        lepton.SFecout = energy / vec_lepton.P();
+        lepton.m2ecout = (m2u + m2v + m2w) / 3;
+      }
+    }
+
+
+    m_log->Debug("Variables obtained from calorimeter bank");
+
+    return lepton;
+  }
+
+
+  double LeptonIDFilter::CalculateScore(LeptonIDVars lepton_vars, concurrent_key_t const key) const
+  {
+    // Assigning variables from lepton_vars for TMVA method
+    std::string weightsfile;
+    switch(lepton_vars.pid) {
+      case 11:
+        weightsfile = o_weightfile_electron->Load(key);
+        break;
+      case -11:
+        weightsfile = o_weightfile_positron->Load(key);
+        break;
+      default:
+        throw std::runtime_error(fmt::format("unknown lepton PDG code {}", lepton_vars.pid));
+    }
+    return readerTMVA->EvaluateMVA(lepton_vars.GetValues(), weightsfile);
+  }
+
+
+  bool LeptonIDFilter::Filter(double score) const
+  {
+    if(score >= o_cut)
+      return true;
+    else
+      return false;
+  }
+
+
+  void LeptonIDFilter::Stop()
+  {
+  }
+
+}