Ft energy correction validator

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

@@ -4,26 +4,26 @@ namespace iguana::clas12 {
   REGISTER_IGUANA_ALGORITHM(FTEnergyCorrection);
 
   void FTEnergyCorrection::Start(hipo::banklist& banks) {
-    b_ft_particle = GetBankIndex(banks, "RECFT::Particle");
+    b_particle = GetBankIndex(banks, "REC::Particle");
     electron_mass = particle::mass.at(particle::electron);
   }
 
   void FTEnergyCorrection::Run(hipo::banklist& banks) const {
-    auto& ftParticleBank = GetBank(banks, b_ft_particle, "RECFT::Particle");
-    ShowBank(ftParticleBank, Logger::Header("INPUT FT PARTICLES"));
-    for(auto const& row : ftParticleBank.getRowList()) {
-      if(ftParticleBank.getInt("pid", row) == particle::PDG::electron) {
-        auto px = ftParticleBank.getFloat("px", row);
-        auto py = ftParticleBank.getFloat("py", row);
-        auto pz = ftParticleBank.getFloat("pz", row);
+    auto& ParticleBank = GetBank(banks, b_particle, "REC::Particle");
+    ShowBank(ParticleBank, Logger::Header("INPUT PARTICLES"));
+    for(auto const& row : ParticleBank.getRowList()) {
+      if(ParticleBank.getInt("pid", row) == particle::PDG::electron) {
+        auto px = ParticleBank.getFloat("px", row);
+        auto py = ParticleBank.getFloat("py", row);
+        auto pz = ParticleBank.getFloat("pz", row);
         auto E = std::hypot(std::hypot(px, py, pz), electron_mass);
         auto [px_new, py_new, pz_new, E_new] = Transform(px, py, pz, E);
-        ftParticleBank.putFloat("px", row, px_new);
-        ftParticleBank.putFloat("py", row, py_new);
-        ftParticleBank.putFloat("pz", row, pz_new);
+        ParticleBank.putFloat("px", row, px_new);
+        ParticleBank.putFloat("py", row, py_new);
+        ParticleBank.putFloat("pz", row, pz_new);
       }
     }
-    ShowBank(ftParticleBank, Logger::Header("OUTPUT FT PARTICLES"));
+    ShowBank(ParticleBank, Logger::Header("OUTPUT PARTICLES"));
   }
 
   vector4_t  FTEnergyCorrection::Transform(

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

@@ -44,7 +44,7 @@ namespace iguana::clas12 {
 
     private:
 
-      hipo::banklist::size_type b_ft_particle;
+      hipo::banklist::size_type b_particle;
       double electron_mass;
 
   };

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

@@ -0,0 +1,160 @@
+#include "Validator.h"
+#include <TProfile.h>
+#include <TStyle.h>
+
+//Edit  m_E_max, -m_deltaE_max and m_deltaE_max.
+//See what the plotting does.
+namespace iguana::clas12 {
+
+  REGISTER_IGUANA_VALIDATOR(FTEnergyCorrectionValidator);
+
+  void FTEnergyCorrectionValidator::Start(hipo::banklist& banks)
+  {
+    // define the algorithm sequence
+    m_algo_seq = std::make_unique<AlgorithmSequence>();
+    m_algo_seq->Add("clas12::FT_Energy_Corr");
+    m_algo_seq->SetOption<std::vector<int>>("clas12::EventBuilderFilter", "pids", u_pdg_list);
+    m_algo_seq->Start(banks);
+
+    // get bank indices
+    b_particle = GetBankIndex(banks, "REC::Particle");
+    electron_mass = particle::mass.at(particle::electron);
+
+    // set an output file
+    auto output_dir = GetOutputDirectory();
+    if(output_dir) {
+      m_output_file_basename = output_dir.value() + "/energy_corrections";
+      m_output_file          = new TFile(m_output_file_basename + ".root", "RECREATE");
+
+    }
+    }
+
+
+  void FTEnergyCorrectionValidator::Run(hipo::banklist& banks) const
+  {
+    // get the momenta before
+    auto& particle_bank = GetBank(banks, b_particle, "REC::Particle");
+    double px_el, py_el, pz_el, E_el;
+    double px_pim, py_pim, pz_pim;
+    double px_pip, py_pip, pz_pip;
+    double px_pr, py_pr, pz_pr;
+    double DeltaE;
+    TLorentzVector beam(0,0,10.6,10.6);
+    TLorentzVector target(0,0,0,0.938272);
+    TLorentzVector miss_el, pion_plus, pion_minus, proton, electron;
+
+
+    for(auto const& row : particle_bank.getRowList()){
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::electron){
+
+            px_el = particle_bank.getFloat("px", row);
+            py_el = particle_bank.getFloat("py", row);
+            pz_el = particle_bank.getFloat("pz", row);
+            electron.SetXYZM(px_el, py_el, pz_el, electron_mass);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::pi_minus){
+
+            px_pim = particle_bank.getFloat("px", row);
+            py_pim = particle_bank.getFloat("py", row);
+            pz_pim = particle_bank.getFloat("pz", row);
+            pion_minus.SetXYZM(px_pim, py_pim, pz_pim, 0.139600);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::pi_plus){
+
+            px_pip = particle_bank.getFloat("px", row);
+            py_pip = particle_bank.getFloat("py", row);
+            pz_pip = particle_bank.getFloat("pz", row);
+            pion_plus.SetXYZM(px_pip, py_pip, pz_pip, 0.139600);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::proton){
+
+            px_pr = particle_bank.getFloat("px", row);
+            py_pr = particle_bank.getFloat("py", row);
+            pz_pr = particle_bank.getFloat("pz", row);
+            proton.SetXYZM(px_pr, py_pr, pz_pr, 0.938272);
+          }
+          miss_el =  beam + target - pion_minus - pion_plus - proton;
+          DeltaE = miss_el.E() - electron.E();
+          h_beforecorr->Fill(electron.E(), DeltaE);
+    }
+    // run the energy corrections
+    m_algo_seq->Run(banks);
+
+    for(auto const& row : particle_bank.getRowList()){
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::electron){
+
+            px_el = particle_bank.getFloat("px", row);
+            py_el = particle_bank.getFloat("py", row);
+            pz_el = particle_bank.getFloat("pz", row);
+            electron.SetXYZM(px_el, py_el, pz_el, electron_mass);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::pi_minus){
+
+            px_pim = particle_bank.getFloat("px", row);
+            py_pim = particle_bank.getFloat("py", row);
+            pz_pim = particle_bank.getFloat("pz", row);
+            pion_minus.SetXYZM(px_pim, py_pim, pz_pim, 0.139600);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::pi_plus){
+
+            px_pip = particle_bank.getFloat("px", row);
+            py_pip = particle_bank.getFloat("py", row);
+            pz_pip = particle_bank.getFloat("pz", row);
+            pion_plus.SetXYZM(px_pip, py_pip, pz_pip, 0.139600);
+          }
+
+          if(particle_bank.getInt("pid", row) == particle::PDG::proton){
+
+            px_pr = particle_bank.getFloat("px", row);
+            py_pr = particle_bank.getFloat("py", row);
+            pz_pr = particle_bank.getFloat("pz", row);
+            proton.SetXYZM(px_pr, py_pr, pz_pr, 0.938272);
+          }
+          miss_el =  beam + target - pion_minus - pion_plus - proton;
+          DeltaE = miss_el.E() - electron.E();
+          h_aftercorr->Fill(electron.E(), DeltaE);
+    }
+
+  }
+  void FTEnergyCorrectionValidator::Stop()
+  {
+    if(GetOutputDirectory()) {
+
+          int n_rows = 1;
+          int n_cols = 2;
+          auto canv = new TCanvas("c","c",n_cols*800,n_rows*800);
+          canv->Divide(n_cols,n_rows);
+          gStyle->SetOptStat(0);
+          for(int pad_num = 0; pad_num < 2; pad_num++){
+              auto pad = canv->GetPad(pad_num+1);
+              pad->cd();
+              pad->SetGrid(1, 1);
+              pad->SetLogz();
+              pad->SetLeftMargin(0.12);
+              pad->SetRightMargin(0.12);
+              pad->SetBottomMargin(0.12);
+
+              switch(pad_num){
+                  case 0:
+                      h_beforecorr->Draw("colz");
+                      break;
+                  case 1:
+                      h_aftercorr->Draw("colz");
+                      break;
+              }
+          }
+      canv->SaveAs(m_output_file_basename + "_plot.png");
+      m_output_file->Write();
+      m_log->Info("Wrote output file {}", m_output_file->GetName());
+      m_output_file->Close();
+    }
+  }
+
+}

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

@@ -0,0 +1,57 @@
+#pragma once
+
+#include "iguana/algorithms/Validator.h"
+#include "iguana/algorithms/TypeDefs.h"
+#include <TH2D.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <Math/Vector4D.h>
+
+namespace iguana::clas12 {
+
+  /// @brief_algo Forward Tagger energy correction
+  ///
+  /// @begin_doc_algo{Transformer}
+  /// @input_banks{RECFT::Particle}
+  /// @output_banks{RECFT::Particle}
+  /// @end_doc
+  class FTEnergyCorrectionValidator : public Validator {
+
+    DEFINE_IGUANA_VALIDATOR(FTEnergyCorrectionValidator, clas12::FTEnergyCorrectionValidator)
+
+    public:
+
+      void Start(hipo::banklist& banks) override;
+      void Run(hipo::banklist& banks) const override;
+      void Stop() override;
+
+      /// @action_function{scalar transformer}
+      /// Transformation function that returns 4-vector of electron with corrected energy for the Forward Tagger.
+      /// Currently only validated for Fall 2018 outbending data.
+      /// @param px @f$p_x@f$
+      /// @param py @f$p_y@f$
+      /// @param pz @f$p_z@f$
+      /// @param E @f$E@f$
+      /// @returns an electron 4-vector with the corrected energy for the Forward Tagger.
+
+      vector4_t Transform(
+          vector_element_t px,
+          vector_element_t py,
+          vector_element_t pz,
+          vector_element_t E) const;
+
+    private:
+
+      // How can I change this FillHistograms function so it works for TH2D histograms? I'm probably being stupid.
+      hipo::banklist::size_type b_particle;
+      double electron_mass;
+      std::vector<int> const u_pdg_list = {particle::PDG::electron, particle::PDG::pi_plus, particle::PDG::pi_minus, particle::PDG::proton};
+      TString m_output_file_basename;
+      TFile* m_output_file;
+
+      TH2D* h_beforecorr;
+      TH2D* h_aftercorr;
+  };
+
+}

src/iguana/algorithms/meson.build

@@ -69,7 +69,6 @@ algo_dict += [
     'name': 'clas12::FTEnergyCorrection',
     'has_config': false,
     'has_c_bindings': false,
-    'has_validator': false,
     'test_args': { 'banks': ['RECFT::Particle'] },
   },
   {