diff --git a/examples/config/my_combined_config_file.yaml b/examples/config/my_combined_config_file.yaml index db7351d1..d26e29bc 100644 --- a/examples/config/my_combined_config_file.yaml +++ b/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: diff --git a/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml b/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml index 1ef2b639..9082fbdb 100644 --- a/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml +++ b/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 ] diff --git a/examples/config/my_z_vertex_cuts.yaml b/examples/config/my_z_vertex_cuts.yaml index 67b0ccaa..a90185af 100644 --- a/examples/config/my_z_vertex_cuts.yaml +++ b/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 ] diff --git a/examples/iguana_ex_fortran_01_action_functions.f b/examples/iguana_ex_fortran_01_action_functions.f index ae4e4164..94d78173 100644 --- a/examples/iguana_ex_fortran_01_action_functions.f +++ b/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 diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc b/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc index fe3aa816..9a872062 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc @@ -11,10 +11,15 @@ namespace iguana::clas12 { ParseYAMLConfig(); o_runnum = GetCachedOption("runnum").value_or(0); // FIXME: should be set form RUN::config o_zcuts = GetOptionVector("cuts", {GetConfig()->InRange("runs", o_runnum), "cuts"}); + o_pids = GetOptionSet("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,8 +37,9 @@ 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; }); @@ -41,9 +47,14 @@ namespace iguana::clas12 { 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(o_pids.find(pid) != o_pids.end()) { + return zvertex > GetZcutLower() && zvertex < GetZcutUpper(); + } else { + return true; + } } int ZVertexFilter::GetRunNum() const diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.h b/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.h index 2a285852..544339fa 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.h +++ b/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 o_zcuts; + + /// pids to apply ZVertexFilter to + std::set o_pids; }; } diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Bindings.cc b/src/iguana/algorithms/clas12/ZVertexFilter/Bindings.cc index 4ff4791e..d448062e 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Bindings.cc +++ b/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(iguana_get_algo_(algo_idx))->Filter(*vz); + *out = *out && dynamic_cast(iguana_get_algo_(algo_idx))->Filter(*vz, *pid); } } } diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml b/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml index ed4cda56..ee02020d 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml @@ -3,12 +3,15 @@ clas12::ZVertexFilter: # default cuts - default: + pids: [ 11 ] cuts: [ -20.0, 20.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 ] diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Validator.cc b/src/iguana/algorithms/clas12/ZVertexFilter/Validator.cc new file mode 100644 index 00000000..a32f74af --- /dev/null +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Validator.cc @@ -0,0 +1,109 @@ +#include "Validator.h" + +#include +#include +#include + +namespace iguana::clas12 { + + REGISTER_IGUANA_VALIDATOR(ZVertexFilterValidator); + + void ZVertexFilterValidator::Start(hipo::banklist& banks) + { + // define the algorithm sequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("clas12::ZVertexFilter"); + m_algo_seq->SetOption>("clas12::ZVertexFilter", "pids", u_pdgtocut_list); + m_algo_seq->SetOption>("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 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 "< 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(); + } + } + +} diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Validator.h b/src/iguana/algorithms/clas12/ZVertexFilter/Validator.h new file mode 100644 index 00000000..f62c36b7 --- /dev/null +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Validator.h @@ -0,0 +1,48 @@ +#pragma once + +#include "iguana/algorithms/TypeDefs.h" +#include "iguana/algorithms/Validator.h" + +#include +#include + +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 const u_pdg_list = { + particle::PDG::electron, + particle::PDG::pi_plus, + particle::PDG::pi_minus, + particle::PDG::proton, + particle::PDG::neutron}; + + std::vector const u_pdgtocut_list = { + particle::PDG::electron}; + + std::vector const u_cuts_list = { + -5, + 5}; + + TString m_output_file_basename; + TFile* m_output_file; + mutable std::unordered_map> u_zvertexplots; + }; + +} diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 6628d396..7ae7415f 100644 --- a/src/iguana/algorithms/meson.build +++ b/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'] }, }, {