ZVertexFilter: Added option for pids and validator

examples/config/my_combined_config_file.yaml

@@ -23,14 +23,17 @@ clas12::ZVertexFilter:
 
   # default cuts
   - default:
+    pids: [ 11 ]
     cuts: [ -33.0, 11.0 ]
 
   #RG-A fall2018 inbending
   - runs: [ 4760, 5419 ]
+    pids: [ 11 ]
     cuts: [ -13.0, 12.0 ]
 
    #RG-A fall2018 outbending
   - runs: [ 5420, 5674 ]
+    pids: [ 11 ]
     cuts: [ -18.0, 10.0 ]
 
 another::Algorithm:

examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml

@@ -3,12 +3,15 @@ clas12::ZVertexFilter:
 
   # default cuts
   - default:
+    pids: [ 11 ]
     cuts: [ -15.0, 15.0 ]
 
   # RG-A fall2018 inbending
   - runs: [ 4760, 5419 ]
+    pids: [ 11 ]
     cuts: [ -5.0, 3.0 ]
 
   # RG-A fall2018 outbending
   - runs: [ 5420, 5674 ]
+    pids: [ 11 ]
     cuts: [ -8.0, 7.0 ]

examples/config/my_z_vertex_cuts.yaml

@@ -3,12 +3,15 @@ clas12::ZVertexFilter:
 
   # default cuts
   - default:
+    pids: [ 11 ]
     cuts: [ -1.5, 1.3 ]
 
   # RG-A fall2018 inbending
   - runs: [ 4760, 5419 ]
+    pids: [ 11 ]
     cuts: [ -0.5, 0.5 ]
 
   # RG-A fall2018 outbending
   - runs: [ 5420, 5674 ]
+    pids: [ 11 ]
     cuts: [ -0.8, 0.7 ]

examples/iguana_ex_fortran_01_action_functions.f

@@ -193,7 +193,7 @@ program iguana_ex_fortran_01_action_functions
           call iguana_clas12_eventbuilderfilter_filter(
      &      algo_eb_filter, pid(i), accept(i))
           call iguana_clas12_zvertexfilter_filter(
-     &      algo_vz_filter, vz(i), accept(i))
+     &      algo_vz_filter, vz(i), pid(i), accept(i))
           print *, '  i = ', i, '  pid = ', pid(i), ' vz = ', vz(i),
      &      '  =>  accept = ', accept(i)
         enddo

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

@@ -11,10 +11,15 @@ namespace iguana::clas12 {
     ParseYAMLConfig();
     o_runnum = GetCachedOption<int>("runnum").value_or(0); // FIXME: should be set form RUN::config
     o_zcuts  = GetOptionVector<double>("cuts", {GetConfig()->InRange("runs", o_runnum), "cuts"});
+    o_pids  = GetOptionVector<int>("pids", {GetConfig()->InRange("runs", o_runnum), "pids"});
     if(o_zcuts.size() != 2) {
       m_log->Error("configuration option 'cuts' must be an array of size 2, but it is {}", PrintOptionValue("cuts"));
       throw std::runtime_error("bad configuration");
     }
+    if(o_pids.size() < 1) {
+      m_log->Error("configuration option 'pids' must have at least one value");
+      throw std::runtime_error("bad configuration");
+    }
 
     // get expected bank indices
     b_particle = GetBankIndex(banks, "REC::Particle");
@@ -32,18 +37,24 @@ namespace iguana::clas12 {
     // filter the input bank for requested PDG code(s)
     particleBank.getMutableRowList().filter([this](auto bank, auto row) {
         auto zvertex = bank.getFloat("vz", row);
-        auto accept  = Filter(zvertex);
-        m_log->Debug("input vz {} -- accept = {}", zvertex, accept);
+        auto pid = bank.getInt("pid", row);
+        auto accept  = Filter(zvertex,pid);
+        m_log->Debug("input vz {} pid {} -- accept = {}", zvertex, pid, accept);
         return accept ? 1 : 0;
         });
 
     // dump the modified bank
     ShowBank(particleBank, Logger::Header("OUTPUT PARTICLES"));
   }
 
-  bool ZVertexFilter::Filter(double const zvertex) const
+  bool ZVertexFilter::Filter(double const zvertex, int pid) const
   {
-    return zvertex > GetZcutLower() && zvertex < GetZcutUpper();
+    //only apply filter if particle pid is in list of pids
+    if(std::find(o_pids.begin(), o_pids.end(), pid) != o_pids.end()) {
+      return zvertex > GetZcutLower() && zvertex < GetZcutUpper();
+    } else {
+      return true;
+    }
   }
 
   int ZVertexFilter::GetRunNum() const

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

@@ -25,10 +25,11 @@ namespace iguana::clas12 {
       void Run(hipo::banklist& banks) const override;
       void Stop() override;
 
-      /// @action_function{scalar filter} checks if the Z Vertex is within specified bounds
+      /// @action_function{scalar filter} checks if the Z Vertex is within specified bounds if pid is one for which the filter should be applied to
       /// @param zvertex the particle Z Vertex to check
+      /// @param pid the particle pid
       /// @returns `true` if `zvertex` is within specified bounds
-      bool Filter(double const zvertex) const;
+      bool Filter(double const zvertex, int pid) const;
 
       /// @returns the current run number
       int GetRunNum() const;
@@ -46,6 +47,9 @@ namespace iguana::clas12 {
 
       /// Z-vertex cut
       std::vector<double> o_zcuts;
+
+      /// pids to apply ZVertexFilter to
+      std::vector<int> o_pids;
   };
 
 }

src/iguana/algorithms/clas12/ZVertexFilter/Bindings.cc

@@ -7,10 +7,11 @@ namespace iguana::bindings {
   /// @see `iguana::clas12::ZVertexFilter::Filter`
   /// @param [in] algo_idx the algorithm index
   /// @param [in] vz
+  /// @param [in] pid
   /// @param [in,out] out the return value
-  void iguana_clas12_zvertexfilter_filter_(algo_idx_t* algo_idx, float* vz, bool* out)
+  void iguana_clas12_zvertexfilter_filter_(algo_idx_t* algo_idx, float* vz, int* pid, bool* out)
   {
-    *out = *out && dynamic_cast<clas12::ZVertexFilter*>(iguana_get_algo_(algo_idx))->Filter(*vz);
+    *out = *out && dynamic_cast<clas12::ZVertexFilter*>(iguana_get_algo_(algo_idx))->Filter(*vz, *pid);
   }
   }
 }

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

@@ -3,12 +3,15 @@ clas12::ZVertexFilter:
 
   # default cuts
   - default:
-    cuts: [ -20.0, 20.0 ]
+    pids: [ 11 ]
+    cuts: [ -20, 20 ]
 
   # RG-A fall2018 inbending
   - runs: [ 4760, 5419 ]
+    pids: [ 11 ]
     cuts: [ -13.0, 12.0 ]
 
   # RG-A fall2018 outbending
   - runs: [ 5420, 5674 ]
+    pids: [ 11 ]
     cuts: [ -18.0, 10.0 ]

src/iguana/algorithms/clas12/ZVertexFilter/Validator.cc

@@ -0,0 +1,109 @@
+#include "Validator.h"
+
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+
+namespace iguana::clas12 {
+
+  REGISTER_IGUANA_VALIDATOR(ZVertexFilterValidator);
+
+  void ZVertexFilterValidator::Start(hipo::banklist& banks)
+  {
+    // define the algorithm sequence
+    m_algo_seq = std::make_unique<AlgorithmSequence>();
+    m_algo_seq->Add("clas12::ZVertexFilter");
+    m_algo_seq->SetOption<std::vector<int>>("clas12::ZVertexFilter", "pids", u_pdgtocut_list);
+    m_algo_seq->SetOption<std::vector<double>>("clas12::ZVertexFilter", "cuts", u_cuts_list);
+    m_algo_seq->Start(banks);
+
+    // get bank indices
+    b_particle = GetBankIndex(banks, "REC::Particle");
+
+    // set an output file
+    auto output_dir = GetOutputDirectory();
+    if(output_dir) {
+      m_output_file_basename = output_dir.value() + "/zvertex_filter";
+      m_output_file          = new TFile(m_output_file_basename + ".root", "RECREATE");
+    }
+
+    // define plots
+    gStyle->SetOptStat(0);
+    for(auto const& pdg : u_pdg_list) {
+      std::vector<TH1D*> zvertexplots;
+      TString particle_name  = particle::name.at(particle::PDG(pdg));
+      TString particle_title = particle::title.at(particle::PDG(pdg));
+      for(int i = 0; i < 2; i++) {
+        TString beforeafter_name  = "before";
+        if (i == 1){ beforeafter_name = "after";}
+
+        zvertexplots.push_back(new TH1D(
+            "zvertexplots_" + particle_name + "_" + beforeafter_name,
+            particle_title + " Z Vertex ; Z Vertex [cm]",
+            200, -40, 40));
+
+        //std::cout<<"Adding plots for "<<pdg<<" "<<beforeafter_name<<std::endl;
+      }
+      u_zvertexplots.insert({pdg, zvertexplots});
+    }
+  }
+
+
+  void ZVertexFilterValidator::Run(hipo::banklist& banks) const
+  {
+    auto& particle_bank = GetBank(banks, b_particle, "REC::Particle");
+
+    // lock the mutex, so we can mutate plots
+    std::scoped_lock<std::mutex> lock(m_mutex);
+
+    // fill the plots before
+    for(auto const& row : particle_bank.getRowList()) {
+      double vz = particle_bank.getFloat("vz", row);
+      int pdg = particle_bank.getInt("pid", row);
+      auto it = u_zvertexplots.find(pdg);
+      //check if pdg is amongs those that we want to plot
+      if (it != u_zvertexplots.end()) {
+        u_zvertexplots.at(pdg).at(0)->Fill(vz);
+      }
+    }
+
+    // run the momentum corrections
+    m_algo_seq->Run(banks);
+
+    // fill the plots after
+    for(auto const& row : particle_bank.getRowList()) {
+      double vz = particle_bank.getFloat("vz", row);
+      int pdg = particle_bank.getInt("pid", row);
+      auto it = u_zvertexplots.find(pdg);
+      //check if pdg is amongs those that we want to plot
+      if (it != u_zvertexplots.end()) {
+        u_zvertexplots.at(pdg).at(1)->Fill(vz);
+      }
+    }
+  }
+
+  void ZVertexFilterValidator::Stop()
+  {
+    if(GetOutputDirectory()) {
+      for(auto const& [pdg, plots] : u_zvertexplots) {
+        TString canv_name = Form("canv%d", pdg);
+        auto canv         = new TCanvas(canv_name, canv_name, 800, 600);
+
+        plots.at(0)->SetLineColor(kBlue);
+        plots.at(0)->SetLineWidth(2);
+        plots.at(0)->Draw("");
+
+        plots.at(1)->SetLineColor(kRed);
+        plots.at(1)->SetLineWidth(2);
+        plots.at(1)->Draw("same");
+
+        canv->Draw();
+        canv->SaveAs(m_output_file_basename + "_" + std::to_string(pdg) + ".png");
+      }
+      m_output_file->Write();
+      m_log->Info("Wrote output file {}", m_output_file->GetName());
+      m_output_file->Close();
+    }
+  }
+
+}

src/iguana/algorithms/clas12/ZVertexFilter/Validator.h

@@ -0,0 +1,48 @@
+#pragma once
+
+#include "iguana/algorithms/TypeDefs.h"
+#include "iguana/algorithms/Validator.h"
+
+#include <TFile.h>
+#include <TH1.h>
+
+namespace iguana::clas12 {
+
+  /// @brief `iguana::clas12::ZVertexFilter` validator
+  class ZVertexFilterValidator : public Validator
+  {
+
+      DEFINE_IGUANA_VALIDATOR(ZVertexFilterValidator, clas12::ZVertexFilterValidator)
+
+    public:
+
+      void Start(hipo::banklist& banks) override;
+      void Run(hipo::banklist& banks) const override;
+      void Stop() override;
+
+    private:
+
+      hipo::banklist::size_type b_particle;
+
+      // add pdgs not to cut to check we're 
+      // only cutting on right particles
+      std::vector<int> const u_pdg_list = {
+          particle::PDG::electron,
+          particle::PDG::pi_plus,
+          particle::PDG::pi_minus,
+          particle::PDG::proton,
+          particle::PDG::neutron};
+
+      std::vector<int> const u_pdgtocut_list = {
+          particle::PDG::electron};
+
+      std::vector<double> const u_cuts_list = {
+          -5,
+          5};
+
+      TString m_output_file_basename;
+      TFile* m_output_file;
+      mutable std::unordered_map<int, std::vector<TH1D*>> u_zvertexplots;
+  };
+
+}

src/iguana/algorithms/meson.build

@@ -49,7 +49,7 @@ algo_dict += [
   },
   {
     'name': 'clas12::ZVertexFilter',
-    'has_validator': false,
+    'has_validator': true,
     'test_args': { 'banks': ['REC::Particle'] },
   },
   {