diff --git a/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.cc index 17fe4784..cca674da 100644 --- a/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.cc +++ b/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( diff --git a/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.h b/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.h index 0e487bc1..022992cb 100644 --- a/src/iguana/algorithms/clas12/FTEnergyCorrection/Algorithm.h +++ b/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; }; diff --git a/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.cc b/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.cc new file mode 100644 index 00000000..54d3992b --- /dev/null +++ b/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.cc @@ -0,0 +1,160 @@ +#include "Validator.h" +#include +#include + +//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(); + m_algo_seq->Add("clas12::FT_Energy_Corr"); + m_algo_seq->SetOption>("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(); + } + } + +} diff --git a/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.h b/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.h new file mode 100644 index 00000000..6ac388e1 --- /dev/null +++ b/src/iguana/algorithms/clas12/FTEnergyCorrection/Validator.h @@ -0,0 +1,57 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" +#include "iguana/algorithms/TypeDefs.h" +#include +#include +#include +#include +#include + +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 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; + }; + +} diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 7d60563a..75ce4595 100644 --- a/src/iguana/algorithms/meson.build +++ b/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'] }, }, {