LeptonID

meson.build

@@ -22,6 +22,14 @@ project_description = 'Implementation Guardian of Analysis Algorithms'
 # initialize binding languanges
 add_languages('fortran', native: false, required: 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')
@@ -124,6 +132,7 @@ if ROOT_dep.found()
     '-lGenVector',
     '-lROOTDataFrame',
     '-lROOTVecOps',
+    '-lTMVA',
     '-lTreePlayer',
   ]
   ROOT_dep_link_args_for_validators = [
@@ -141,6 +150,7 @@ project_inc       = include_directories('src')
 project_libs      = []
 project_deps      = declare_dependency(dependencies: [ fmt_dep, yamlcpp_dep, hipo_dep ] ) # do NOT include ROOT here
 project_etc       = get_option('sysconfdir') / meson.project_name()
+project_data      = get_option('datadir') / meson.project_name()
 project_test_env  = environment()
 project_pkg_vars  = [
   'bindir='            + '${prefix}' / get_option('bindir'),
@@ -171,7 +181,9 @@ project_test_env.set(
 
 # set preprocessor macros
 add_project_arguments(
-  '-DIGUANA_ETCDIR="' + get_option('prefix') / project_etc + '"',
+  '-DIGUANA_PREFIX="' + get_option('prefix') + '"',
+  '-DIGUANA_ETCDIR="' + project_etc + '"',
+  '-DIGUANA_DATADIR="' + project_data + '"',
   language: ['cpp'],
 )
 

meson/minimum-version.sh

@@ -32,7 +32,7 @@ case $dep in
     [ "$cmd" = "src" ] && echo "ERROR: command '$cmd' is not used for '$dep'" >&2 && exit 1
     ;;
   root|ROOT)
-    result_meson='>=6.28'
+    result_meson='>=6.26'
     [ "$cmd" = "ALA" ] && echo "ERROR: command '$cmd' is not used for '$dep'" >&2 && exit 1
     result_src='https://root.cern/download/root_v6.28.12.source.tar.gz'
     ;;

src/iguana/algorithms/Algorithm.cc

@@ -112,6 +112,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()
   {
     if(!m_yaml_config) {

src/iguana/algorithms/Algorithm.h

@@ -8,6 +8,7 @@
 
 #include "iguana/algorithms/AlgorithmBoilerplate.h"
 #include "iguana/services/YAMLReader.h"
+#include "iguana/services/DataFileReader.h"
 
 namespace iguana {
 
@@ -136,6 +137,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`.
@@ -226,6 +232,9 @@ namespace iguana {
 
       /// YAML reader
       std::unique_ptr<YAMLReader> m_yaml_config;
+
+      /// Data file reader
+      std::unique_ptr<DataFileReader> m_datafile_reader;
   };
 
   //////////////////////////////////////////////////////////////////////////////

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

@@ -0,0 +1,152 @@
+#include "Algorithm.h"
+
+#include <cmath>
+#include <Math/Vector4D.h>
+
+namespace iguana::clas12 {
+
+  REGISTER_IGUANA_ALGORITHM(LeptonIDFilter , "clas12::LeptonIDFilter");
+
+  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()
+  {
+  }
+
+}

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

@@ -0,0 +1,137 @@
+#pragma once
+
+#include "iguana/algorithms/Algorithm.h"
+#include <TMVA/Reader.h>
+
+///Struct to store variables 
+    struct LeptonIDVars {
+        /// @brief momentum
+        double P;
+        /// @brief Theta angle
+        double Theta;
+        /// @brief Phi angle
+        double Phi;
+        /// @brief Sampling fraction on the PCAL
+        double SFpcal;
+        /// @brief Sampling fraction on the ECIN
+        double SFecin;
+        /// @brief Sampling fraction on the ECOUT
+        double SFecout;
+        /// @brief Second-momenta of PCAL
+        double m2pcal;
+        /// @brief Second-momenta of ECIN
+        double m2ecin;
+        /// @brief Second-momenta of ECOUT
+        double m2ecout;
+        /// @brief Score
+        double score;
+    };
+
+namespace iguana::clas12 {
+  ///
+  /// @brief_algo Filter the leptons from the pion contamination using TMVA models
+  ///
+  /// For each lepton, either positron or electron, it takes some variables from `REC::Particle` (P, Theta and Phi) and `REC::Particle` (Sampling fraction and second moments).
+  /// Using those variables, it call the TMVA method using the weight file, and it computes a score. By a pplying a cut to the score we can separate leptons from pions.
+  ///
+  /// @begin_doc_algo{clas12::LeptonIDFilter | Filter}
+  /// @input_banks{REC::Particle,REC::Calorimeter}
+  /// @end_doc
+  ///
+  /// @begin_doc_config
+  /// @config_param{o_pid | int | PID of the particle; -11 for positrons and 11 for electrons}
+  /// @config_param{o_weightfile | std::string | Location of the weight file of the classifier}
+  /// @config_param{o_cut | double | Value of the score to apply the cut. The algorith will keep all particles that have a score grater than ths value}
+  /// @end_doc
+  class LeptonIDFilter : public Algorithm
+  {
+
+      DEFINE_IGUANA_ALGORITHM(LeptonIDFilter, clas12::LeptonIDFilter)
+
+    public:
+
+      void Start(hipo::banklist& banks) override;
+      void Run(hipo::banklist& banks) const override;
+      void Stop() override;
+
+      /// **FindLepton function**: returns the pindex of the lepton
+      /// @param particle_bank the particle bank
+      /// @returns pindex of the lepton, -1 if there is no lepton
+      int FindLepton(hipo::bank const& particle_bank) const;
+
+
+      /// **GetLeptonIDVariables function**: Using the pindex retrieves the necessary variables from banks
+      /// @param plepton pindex of the lepton
+      /// @param particle_bank the particle bank
+      /// @param calorimeter_bank the calorimeter bank
+      /// @returns LeptonIDVars, the variables required for identification
+      LeptonIDVars GetLeptonIDVariables(int const plepton, hipo::bank const& particle_bank, hipo::bank const& calorimeter_bank) const;
+
+
+      /// **CalculateScore function**: Using the LeptonIDVars variables calculate the score
+      /// @param lepton_vars LeptonIDVars variables
+      /// @returns double, the score
+      double CalculateScore(LeptonIDVars lepton_vars) const;
+
+      /// **Filter function**: Returns true if the particle passed the cut 
+      /// @param score the score obtained from the CalculateScore function
+      /// @returns bool, true if score>=cut, false otherwise
+      bool Filter(double score) const;
+
+      //Create TMVA reader
+      TMVA::Reader *readerTMVA = new TMVA::Reader();
+
+      ///Set of variables for the reader
+      ///Momentum
+      Float_t P;
+      ///Theta angle
+      Float_t Theta;
+      ///Phi angle
+      Float_t Phi;
+      ///Sampling fraction on the PCAL
+      Float_t PCAL;
+      ///Sampling fraction on the ECIN
+      Float_t ECIN;
+      ///Sampling fraction on the ECOUT
+      Float_t ECOUT;
+      ///Second-momenta of PCAL
+      Float_t m2PCAL;
+      ///Second-momenta of ECIN
+      Float_t m2ECIN;
+      ///Second-momenta of ECOUT
+      Float_t m2ECOUT;
+
+      /// @brief Add variables to the readerTMVA
+      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);
+
+      readerTMVA->BookMVA( "BDT", o_weightfile_fullpath );
+
+
+    private:
+
+
+      /// `hipo::banklist` 
+      hipo::banklist::size_type b_particle;
+      hipo::banklist::size_type b_calorimeter;
+
+
+
+
+      /// pid of the lepton
+      int o_pid;
+      /// Location of the weight file
+      std::string o_weightfile;
+      std::string o_weightfile_fullpath;
+      /// Value of the cut
+      double o_cut;
+  };
+
+}

src/iguana/algorithms/clas12/LeptonIDFilter/Config.yaml

@@ -0,0 +1,4 @@
+clas12::LeptonIDFilter:
+  pid: -11
+  weightfile: "weights/9_BDT_positrons_S19.weights.xml"
+  cut: 0.0