LeptonID

doc/setup.md

@@ -1,10 +1,11 @@
 # Setup Guide
 
-| **Table of Contents**                           |
-| ---                                             |
-| 🟠 [**Dependencies**](#dependencies)            |
-| 🟠 [**Building and Installing**](#building)     |
-| 🟠 [**Environment Variables** (optional)](#env) |
+| **Table of Contents**                                      |
+| ---                                                        |
+| 🟠 [**Dependencies**](#dependencies)                       |
+| 🟠 [**Building and Installing**](#building)                |
+| 🟠 [**Environment Variables** (optional)](#env)            |
+| 🟠 [**Connect to a C++ Analysis (consuming)**](#consuming) |
 
 <a name="dependencies"></a>
 ## 🟠 Dependencies
@@ -209,5 +210,23 @@ The following environment variables are set or modified; not all of them are nee
 | `LD_LIBRARY_PATH` (Linux) or `DYLD_LIBRARY_PATH` (macOS) | adds paths to dependency and Iguana libraries                                                                                             |
 | `PYTHONPATH`                                             | adds paths to dependency and Iguana Python packages, if Python bindings are installed                                                     |
 | `ROOT_INCLUDE_PATH`                                      | adds paths to dependency and Iguana header files, for usage in ROOT                                                                       |
-| `IGUANA_CONFIG_PATH`                                     | path to iguana algorithm configuration files (`.yaml`); users may override this with their own path; multiple paths may be specified, delimited by colons (`:`), where paths listed first will override paths listed later (similar behavior as `$PATH`); this variable is _only_ necessary if the Iguana installation has been relocated |
-| `IGUANA`                                                 | the path to the Iguana installation prefix, equivalent to `pkg-config iguana --variable prefix`; this is only for consumers that do not use `pkg-config` or the other standard environment variables, however usage of this variable is _discouraged_ since the installation layout may vary |
+| `IGUANA_CONFIG_PATH`                                     | path to iguana algorithm configuration files (`.yaml`); users may override this with their own path; multiple paths may be specified, delimited by colons (`:`), where paths listed first will override paths listed later (similar behavior as `$PATH`) |
+| `IGUANA`                                                 | the path to the Iguana installation prefix |
+
+> [!IMPORTANT]
+> If the Iguana installation is _relocated_, the following environment variables become _necessary_:
+> - `IGUANA`
+> - `IGUANA_CONFIG_PATH`
+
+> [!WARNING]
+> Avoid using the `$IGUANA` environment variable for locating Iguana libraries, headers, _etc._, since the installation
+> tree may vary. Instead, use `pkg-config`, following the [consuming section](#consuming).
+
+<a name="consuming"></a>
+## 🟠 Connect to a C++ Analysis (consuming)
+
+To integrate Iguana with your own C++ analysis code, _i.e._ to "consume" Iguana as a dependency, follow
+one of the examples, depending on your build tool:
+- [CMake](../examples/build_with_cmake)
+- [Meson](../examples/build_with_meson)
+- [Makefile](../examples/build_with_make)

examples/meson.build

@@ -1,5 +1,5 @@
 # install config files
-example_config_files_prefix = project_etc / '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

@@ -8,6 +8,7 @@ project(
     'cpp_std': 'c++17',
     'buildtype': 'release',
     'libdir': 'lib',
+    'datadir': 'share',
     'licensedir': 'share/licenses/iguana',
     'pkgconfig.relocatable': 'true',
     'force_fallback_for': ['rcdb'],
@@ -24,6 +25,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
+
 # meson modules
 pkg = import('pkgconfig')
 fs = import('fs')
@@ -150,6 +159,7 @@ if ROOT_dep.found()
     '-lGenVector',
     '-lROOTDataFrame',
     '-lROOTVecOps',
+    '-lTMVA',
     '-lTreePlayer',
   ]
   # additional ROOT libraries for validators (namely, graphics libraries)
@@ -163,12 +173,13 @@ if ROOT_dep.found()
 endif
 
 # general project vars
-project_inc = include_directories('src')
-project_libs = []
-project_deps = declare_dependency(dependencies: dep_list) # do NOT include ROOT here
-project_etc = get_option('sysconfdir') / meson.project_name()
-project_test_env = environment()
-project_pkg_vars = [
+project_inc       = include_directories('src')
+project_libs      = []
+project_deps      = declare_dependency(dependencies: dep_list) # do NOT include ROOT here
+project_etc_dir   = get_option('sysconfdir') / meson.project_name()
+project_data_dir  = get_option('datadir') / meson.project_name()
+project_test_env  = environment()
+project_pkg_vars  = [
   'bindir='            + '${prefix}' / get_option('bindir'),
   'dep_pkgconfigdirs=' + ':'.join(dep_pkgconfig_dirs),
   'dep_libdirs='       + ':'.join(dep_lib_dirs),
@@ -198,8 +209,10 @@ project_test_env.set(
 
 # set preprocessor macros
 add_project_arguments(
-  '-DIGUANA_ETCDIR="' + get_option('prefix') / project_etc + '"',
-  language: [ 'cpp' ],
+  '-DIGUANA_PREFIX="' + get_option('prefix') + '"',
+  '-DIGUANA_ETCDIR="' + project_etc_dir + '"',
+  '-DIGUANA_DATADIR="' + project_data_dir + '"',
+  language: ['cpp'],
 )
 
 # start chameleon
@@ -257,7 +270,7 @@ if get_option('z_install_envfile')
       'python': get_option('bind_python') ? 'true' : 'false',
       'libdir': get_option('libdir'),
       'root': ROOT_dep.found() ? 'true' : 'false',
-      'etcdir': project_etc,
+      'etcdir': project_etc_dir,
     },
   )
   foreach shell : [ 'csh', 'tcsh' ]

src/iguana/algorithms/Algorithm.cc

@@ -106,6 +106,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

@@ -9,7 +9,8 @@
 
 #include "iguana/algorithms/AlgorithmBoilerplate.h"
 #include "iguana/services/YAMLReader.h"
-#include <iguana/services/GlobalParam.h>
+#include "iguana/services/GlobalParam.h"
+#include "iguana/services/DataFileReader.h"
 
 namespace iguana {
 
@@ -135,6 +136,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);
+
     protected: // methods
 
       /// Parse YAML configuration files. Sets `m_yaml_config`.
@@ -231,7 +237,8 @@ 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,166 @@
+#include "Algorithm.h"
+
+#include <cmath>
+#include <Math/Vector4D.h>
+
+namespace iguana::clas12 {
+
+  REGISTER_IGUANA_ALGORITHM(LeptonIDFilter , "clas12::LeptonIDFilter");
+
+void LeptonIDFilter::initializeTMVA() {
+          readerTMVA = std::make_unique<TMVA::Reader>("V");
+          // Initialize the variables for the TMVA reader
+          readerTMVA->AddVariable( "P", &P );
+          readerTMVA->AddVariable( "Theta", &Theta );
+          readerTMVA->AddVariable( "Phi", &Phi );
+          readerTMVA->AddVariable( "SFPCAL", &PCAL );
+          readerTMVA->AddVariable( "SFECIN", &ECIN );
+          readerTMVA->AddVariable( "SFECOUT", &ECOUT );
+          readerTMVA->AddVariable( "m2PCAL", &m2PCAL );
+          readerTMVA->AddVariable( "m2ECIN", &m2ECIN );
+          readerTMVA->AddVariable( "m2ECOUT", &m2ECOUT );
+      }
+
+  void LeptonIDFilter::Start(hipo::banklist& banks)
+  {
+    //Get configuration
+    ParseYAMLConfig();
+    o_pid        = GetOptionScalar<int>("pid");//Obtain pid from config file (+11/-11)
+    o_weightfile = GetOptionScalar<std::string>("weightfile");//Obtain weightfile from config file
+    o_cut        = GetOptionScalar<double>("cut");
+
+    // load the weights file
+    o_weightfile_fullpath = GetDataFile(o_weightfile);
+    m_log->Debug("Loaded weight file {}", o_weightfile_fullpath);
+
+    //Get Banks that we are going to use
+    b_particle = GetBankIndex(banks, "REC::Particle");
+    b_calorimeter = GetBankIndex(banks, "REC::Calorimeter");
+
+
+  }
+
+
+  void LeptonIDFilter::Run(hipo::banklist& banks) const
+  {
+    auto& particleBank = GetBank(banks, b_particle, "REC::Particle");
+    auto& calorimeterBank = GetBank(banks, b_calorimeter, "REC::Calorimeter");
+
+    ShowBank(particleBank, Logger::Header("INPUT PARTICLES"));
+
+    //
+    particleBank.getMutableRowList().filter([this,&particleBank,&calorimeterBank](auto bank, auto row) {
+        auto lepton_pindex = FindLepton(particleBank);
+        auto lepton_vars=GetLeptonIDVariables(lepton_pindex,particleBank,calorimeterBank);
+        lepton_vars.score=CalculateScore(lepton_vars);
+
+        return Filter(lepton_vars.score);
+        });
+
+    // dump the modified bank
+    ShowBank(particleBank, Logger::Header("OUTPUT PARTICLES"));
+
+  }
+
+
+  int LeptonIDFilter::FindLepton(hipo::bank const& particle_bank) const{  
+    int lepton_pindex= -1;
+    for(int row = 0; row < particle_bank.getRows(); row++) {
+        auto status = particle_bank.getShort("status", row);
+        if(particle_bank.getInt("pid", row) == o_pid && abs(status)>=2000 && abs(status)<4000) {
+          lepton_pindex=row;
+          break;
+        }
+      }
+      if(lepton_pindex >= 0)
+        m_log->Debug("Found lepton: pindex={}", lepton_pindex);
+      else
+        m_log->Debug("Lepton not found");
+      return lepton_pindex;
+  }
+
+  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) const{
+
+      ///Assing variables from lepton_vars for TMVA method
+      P=lepton_vars.P;
+      Theta=lepton_vars.Theta;
+      Phi=lepton_vars.Phi;
+      PCAL=lepton_vars.SFpcal;
+      ECIN=lepton_vars.SFecin;
+      ECOUT=lepton_vars.SFecout;
+      m2PCAL=lepton_vars.m2pcal;
+      m2ECIN=lepton_vars.m2ecin;
+      m2ECOUT=lepton_vars.m2ecout;
+
+      m_log->Debug("Add variables to readerTMVA");
+      auto score=readerTMVA->EvaluateMVA("BDT");  
+
+      return score;
+  }
+
+  bool LeptonIDFilter::Filter(double score) const{
+    if(score>=o_cut)
+      return true;
+    else
+      return false;
+  }
+
+
+
+  void LeptonIDFilter::Stop()
+  {
+  }
+
+}