diff --git a/doc/gen/Doxyfile.in b/doc/gen/Doxyfile.in index fad3d8d9..c923bc87 100644 --- a/doc/gen/Doxyfile.in +++ b/doc/gen/Doxyfile.in @@ -302,6 +302,8 @@ ALIASES += end_doc_example="@}" # ignore things we don't want to document ALIASES += doxygen_off="@cond NO_DOC" ALIASES += doxygen_on="@endcond" +# latex +ALIASES += latex{1}="@f$\1@f$" # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For diff --git a/examples/iguana_ex_fortran_01_action_functions.f b/examples/iguana_ex_fortran_01_action_functions.f index 92a7766f..d3950c6b 100644 --- a/examples/iguana_ex_fortran_01_action_functions.f +++ b/examples/iguana_ex_fortran_01_action_functions.f @@ -51,6 +51,7 @@ program iguana_ex_fortran_01_action_functions logical(c_bool) accept(N_MAX) ! filter real(c_double) qx, qy, qz, qE ! q vector real(c_double) Q2, x, y, W, nu ! inclusive kinematics + real(c_double) beamPz, targetM ! beam and target integer(c_int) key_vz_filter ! key for Z-vertex filter integer(c_int) key_inc_kin ! key for inclusive kinematics @@ -273,7 +274,8 @@ program iguana_ex_fortran_01_action_functions & px(i_ele), py(i_ele), pz(i_ele), & key_inc_kin, & qx, qy, qz, qE, - & Q2, x, y, W, nu) + & Q2, x, y, W, nu, + & beamPz, targetM) print *, '===> inclusive kinematics:' print *, ' q = (', qx, qy, qz, qE, ')' print *, ' Q2 = ', Q2 diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index a15f8838..49331342 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -74,6 +74,15 @@ algo_dict = [ 'algorithm_needs_ROOT': true, 'test_args': {'banks': [ 'REC::Particle', 'RUN::config' ]}, }, + { + 'name': 'physics::DihadronKinematics', + 'algorithm_needs_ROOT': true, + 'has_action_yaml': false, # FIXME: needs a vector action function + 'test_args': { + 'banks': [ 'REC::Particle', 'RUN::config' ], + 'prerequisites': [ 'physics::InclusiveKinematics' ], + }, + }, { 'name': 'clas12::PhotonGBTFilter', 'algorithm_needs_ROOT': true, diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc new file mode 100644 index 00000000..6b8a95b4 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc @@ -0,0 +1,298 @@ +#include "Algorithm.h" +#include "TypeDefs.h" + +#include + +namespace iguana::physics { + + REGISTER_IGUANA_ALGORITHM(DihadronKinematics, "physics::DihadronKinematics"); + + void DihadronKinematics::Start(hipo::banklist& banks) + { + b_particle = GetBankIndex(banks, "REC::Particle"); + b_inc_kin = GetBankIndex(banks, "physics::InclusiveKinematics"); + + // create the output bank + // FIXME: generalize the groupid and itemid + auto result_schema = CreateBank( + banks, + b_result, + GetClassName(), + { + "pindex_a/S", + "pindex_b/S", + "pdg_a/I", + "pdg_b/I", + "Mh/D", + "z/D", + "MX/D", + "xF/D", + "phiH/D", + "phiR/D", + "theta/D" + }, + 0xF000, + 5); + i_pindex_a = result_schema.getEntryOrder("pindex_a"); + i_pindex_b = result_schema.getEntryOrder("pindex_b"); + i_pdg_a = result_schema.getEntryOrder("pdg_a"); + i_pdg_b = result_schema.getEntryOrder("pdg_b"); + i_Mh = result_schema.getEntryOrder("Mh"); + i_z = result_schema.getEntryOrder("z"); + i_MX = result_schema.getEntryOrder("MX"); + i_xF = result_schema.getEntryOrder("xF"); + i_phiH = result_schema.getEntryOrder("phiH"); + i_phiR = result_schema.getEntryOrder("phiR"); + i_theta = result_schema.getEntryOrder("theta"); + + // parse config file + ParseYAMLConfig(); + o_hadron_a_pdgs = GetOptionSet("hadron_a_list"); + o_hadron_b_pdgs = GetOptionSet("hadron_b_list"); + o_phi_r_method = GetOptionScalar("phi_r_method"); + o_theta_method = GetOptionScalar("theta_method"); + + // check phiR method + if(o_phi_r_method == "RT_via_covariant_kT") + m_phi_r_method = e_RT_via_covariant_kT; + else + throw std::runtime_error(fmt::format("unknown phi_r_method: {:?}", o_phi_r_method)); + + // check theta method + if(o_theta_method == "hadron_a") + m_theta_method = e_hadron_a; + else + throw std::runtime_error(fmt::format("unknown theta_method: {:?}", o_theta_method)); + + } + + /////////////////////////////////////////////////////////////////////////////// + + void DihadronKinematics::Run(hipo::banklist& banks) const + { + auto& particle_bank = GetBank(banks, b_particle, "REC::Particle"); + auto& inc_kin_bank = GetBank(banks, b_inc_kin, "physics::InclusiveKinematics"); + auto& result_bank = GetBank(banks, b_result, GetClassName()); + ShowBank(particle_bank, Logger::Header("INPUT PARTICLES")); + + if(particle_bank.getRowList().empty() || inc_kin_bank.getRowList().empty()) { + m_log->Debug("skip this event, since not all required banks have entries"); + return; + } + + // get beam and target momenta + // FIXME: makes some assumptions about the beam; this should be generalized... + ROOT::Math::PxPyPzMVector p_beam( + 0.0, + 0.0, + inc_kin_bank.getDouble("beamPz", 0), + particle::mass.at(particle::electron)); + ROOT::Math::PxPyPzMVector p_target( + 0.0, + 0.0, + 0.0, + inc_kin_bank.getDouble("targetM", 0)); + + // get virtual photon momentum + ROOT::Math::PxPyPzEVector p_q( + inc_kin_bank.getDouble("qx", 0), + inc_kin_bank.getDouble("qy", 0), + inc_kin_bank.getDouble("qz", 0), + inc_kin_bank.getDouble("qE", 0)); + + // get additional inclusive variables + auto W = inc_kin_bank.getDouble("W", 0); + + // boosts + ROOT::Math::Boost boost__qp((p_q + p_target).BoostToCM()); // CoM frame of target and virtual photon + auto p_q__qp = boost__qp(p_q); + + // build list of dihadron rows (pindices) + auto dih_rows = PairHadrons(particle_bank); + + // loop over dihadrons + result_bank.setRows(dih_rows.size()); + int dih_row = 0; + for(const auto& [row_a, row_b] : dih_rows) { + + // get hadron momenta + Hadron had_a{.row = row_a}; + Hadron had_b{.row = row_b}; + for(auto& had : {&had_a, &had_b}) { + had->pdg = particle_bank.getInt("pid", had->row); + had->p = ROOT::Math::PxPyPzMVector( + particle_bank.getFloat("px", had->row), + particle_bank.getFloat("py", had->row), + particle_bank.getFloat("pz", had->row), + particle::mass.at(static_cast(had->pdg))); + } + + // calculate dihadron momenta and boosts + auto p_Ph = had_a.p + had_b.p; + auto p_Ph__qp = boost__qp(p_Ph); + ROOT::Math::Boost boost__dih(p_Ph.BoostToCM()); // CoM frame of dihadron + + // calculate z + double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); + + // calculate Mh + double Mh = p_Ph.M(); + + // calculate MX + double MX = (p_target + p_q - p_Ph).M(); + + // calculate xF + double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); + + // calculate phiH + double phiH = PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + p_Ph.Vect()).value_or(UNDEF); + + // calculate PhiR + double phiR = UNDEF; + switch(m_phi_r_method) { + case e_RT_via_covariant_kT: + { + for(auto& had : {&had_a, &had_b}) { + had->z = p_target.Dot(had->p) / p_target.Dot(p_q); + had->p_perp = RejectVector(had->p.Vect(), p_q.Vect()); + } + if(had_a.p_perp.has_value() && had_b.p_perp.has_value()) { + auto RT = (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()) / (had_a.z + had_b.z); + phiR = PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + RT).value_or(UNDEF); + } + break; + } + } + + // calculate theta + double theta = UNDEF; + switch(m_theta_method) { + case e_hadron_a: + { + theta = VectorAngle( + boost__dih(had_a.p).Vect(), + p_Ph.Vect()).value_or(UNDEF); + break; + } + } + + result_bank.putShort(i_pindex_a, dih_row, had_a.row); + result_bank.putShort(i_pindex_b, dih_row, had_b.row); + result_bank.putInt(i_pdg_a, dih_row, had_a.pdg); + result_bank.putInt(i_pdg_b, dih_row, had_b.pdg); + result_bank.putDouble(i_Mh, dih_row, Mh); + result_bank.putDouble(i_z, dih_row, z); + result_bank.putDouble(i_MX, dih_row, MX); + result_bank.putDouble(i_xF, dih_row, xF); + result_bank.putDouble(i_phiH, dih_row, phiH); + result_bank.putDouble(i_phiR, dih_row, phiR); + result_bank.putDouble(i_theta, dih_row, theta); + + dih_row++; + } + + ShowBank(result_bank, Logger::Header("CREATED BANK")); + } + + /////////////////////////////////////////////////////////////////////////////// + + std::vector> DihadronKinematics::PairHadrons(hipo::bank const& particle_bank) const { + std::vector> result; + // loop over REC::Particle rows, for hadron A + for(auto const& row_a : particle_bank.getRowList()) { + // check PDG is in the hadron-A list + if(auto pdg_a{particle_bank.getInt("pid", row_a)}; o_hadron_a_pdgs.find(pdg_a) != o_hadron_a_pdgs.end()) { + // loop over REC::Particle rows, for hadron B + for(auto const& row_b : particle_bank.getRowList()) { + // don't pair a particle with itself + if(row_a == row_b) + continue; + // check PDG is in the hadron-B list + if(auto pdg_b{particle_bank.getInt("pid", row_b)}; o_hadron_b_pdgs.find(pdg_b) != o_hadron_b_pdgs.end()) { + // if the PDGs of hadrons A and B are the same, don't double count + if(pdg_a == pdg_b && row_b < row_a) + continue; + // we have a unique dihadron, add it to the list + result.push_back({row_a, row_b}); + } + } + } + } + // trace logging + if(m_log->GetLevel() <= Logger::Level::trace) { + if(result.empty()) + m_log->Trace("=> no dihadrons in this event"); + else + m_log->Trace("=> number of dihadrons found: {}", result.size()); + } + return result; + } + + /////////////////////////////////////////////////////////////////////////////// + + std::optional DihadronKinematics::PlaneAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b, + ROOT::Math::XYZVector const v_c, + ROOT::Math::XYZVector const v_d) + { + auto cross_ab = v_a.Cross(v_b); // A x B + auto cross_cd = v_c.Cross(v_d); // C x D + + auto sgn = cross_ab.Dot(v_d); // (A x B) . D + if(!(std::abs(sgn) > 0)) + return std::nullopt; + sgn /= std::abs(sgn); // sign of (A x B) . D + + auto numer = cross_ab.Dot(cross_cd); // (A x B) . (C x D) + auto denom = cross_ab.R() * cross_cd.R(); // |A x B| * |C x D| + if(!(std::abs(denom) > 0)) + return std::nullopt; + return sgn * std::acos(numer / denom); + } + + std::optional DihadronKinematics::ProjectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + auto denom = v_b.Dot(v_b); + if(!(std::abs(denom) > 0)) + return std::nullopt; + return v_b * ( v_a.Dot(v_b) / denom ); + } + + std::optional DihadronKinematics::RejectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + auto v_c = ProjectVector(v_a, v_b); + if(v_c.has_value()) + return v_a - v_c.value(); + return std::nullopt; + } + + std::optional DihadronKinematics::VectorAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + double m = v_a.R() * v_b.R(); + if(m > 0) + return std::acos(v_a.Dot(v_b) / m); + return std::nullopt; + } + + /////////////////////////////////////////////////////////////////////////////// + + void DihadronKinematics::Stop() + { + } + +} diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h new file mode 100644 index 00000000..2a761b2e --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h @@ -0,0 +1,163 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" +#include +#include + +namespace iguana::physics { + + /// Set of dihadron kinematics variables + struct DihadronKinematicsVars { + /// @brief `REC::Particle` row (`pindex`) of hadron A + int pindex_a; + /// @brief `REC::Particle` row (`pindex`) of hadron B + int pindex_b; + /// @brief PDG code of hadron A + int pdg_a; + /// @brief PDG code of hadron B + int pdg_b; + /// @brief @latex{M_h}: Invariant mass of the dihadron + double Mh; + /// @brief @latex{z}: Fraction of energy of fragmenting parton carried by the dihadron + double z; + /// @brief @latex{M_X(ehhX)}: Missing mass of the dihadron + double MX; + /// @brief @latex{x_F}: Feynman-x of the dihadron + double xF; + /// @brief @latex{\phi_h}: @latex{q}-azimuthal angle between the lepton-scattering plane and the @latex{\vec{q}\times\vec{P}_h} plane; + /// if the value is `DihadronKinematics::UNDEF`, the calculation failed + double phiH; + /// @brief @latex{\phi_R}: @latex{q}-azimuthal angle between the lepton-scattering plane and dihadron plane; + /// if the value is `DihadronKinematics::UNDEF`, the calculation failed + double phiR; + /// @brief @latex{\theta}: the "decay" angle of hadron A in the dihadron rest frame, with respect; + /// to the dihadron momentum direction + double theta; + }; + + /// @brief_algo Calculate semi-inclusive dihadron kinematic quantities defined in `iguana::physics::DihadronKinematicsVars` + /// + /// @begin_doc_algo{physics::DihadronKinematics | Creator} + /// @input_banks{REC::Particle, physics::InclusiveKinematics} + /// @output_banks{%physics::DihadronKinematics} + /// @end_doc + /// + /// @begin_doc_config + /// @config_param{hadron_a_list | list[int] | list of "hadron A" PDGs} + /// @config_param{hadron_b_list | list[int] | list of "hadron B" PDGs} + /// @config_param{phi_r_method | string | method used to calculate @latex{\phi_R} (see section "phiR calculation methods" below)} + /// @config_param{theta_method | string | method used to calculate @latex{\theta} (see section "theta calculation methods" below)} + /// @end_doc + /// + /// Dihadron PDGs will be formed from pairs from `hadron_a_list` and `hadron_b_list`. For example, + /// if you define: + /// ```yaml + /// hadron_a_list: [ 211 ] + /// hadron_b_list: [ -211, 2212 ] + /// ``` + /// then the algorithm will calculate kinematics for @latex{\pi^+\pi^-} and @latex{\pi^+p} dihadrons; hadron A + /// is the @latex{\pi^+} for both of these, whereas hadron B is the @latex{\pi^-} for the former and the proton + /// for the latter. + /// + /// @par phiR calculation methods + /// - `"RT_via_covariant_kT"`: use @latex{R_T} computed via covariant @latex{k_T} formula + /// + /// @par theta calculation methods + /// - `"hadron_a"`: use hadron A's "decay angle" in the dihadron rest frame + class DihadronKinematics : public Algorithm + { + + DEFINE_IGUANA_ALGORITHM(DihadronKinematics, physics::DihadronKinematics) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + /// a value used when some calculation fails + double const UNDEF{-10000}; + + /// @brief form dihadrons by pairing hadrons + /// @param particle_bank the particle bank (`REC::Particle`) + /// @returns a list of pairs of hadron rows + std::vector> PairHadrons(hipo::bank const& particle_bank) const; + + /// @brief calculate the angle between two planes + /// + /// The two planes are transverse to @latex{\vec{v}_a\times\vec{v}_b} and @latex{\vec{v}_c\times\vec{v}_d} + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @param v_c vector @latex{\vec{v}_c} + /// @param v_d vector @latex{\vec{v}_d} + /// @returns the angle between the planes, in radians, if the calculation is successful + static std::optional PlaneAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b, + ROOT::Math::XYZVector const v_c, + ROOT::Math::XYZVector const v_d); + + /// @brief projection of one vector onto another + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the vector @latex{\vec{v}_a} projected onto vector @latex{\vec{v}_b}, if the calculation is successful + static std::optional ProjectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + + /// @brief projection of one vector onto the plane transverse to another vector + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the vector @latex{\vec{v}_a} projected onto the plane transverse to @latex{\vec{v}_b}, if the calculation is successful + static std::optional RejectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + + /// @brief calculate the angle between two vectors + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the angle between @latex{\vec{v}_a} and @latex{\vec{v}_b}, if the calculation is successful + static std::optional VectorAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + + private: + + // banklist indices + hipo::banklist::size_type b_particle; + hipo::banklist::size_type b_inc_kin; + hipo::banklist::size_type b_result; + + // `b_result` bank item indices + int i_pindex_a; + int i_pindex_b; + int i_pdg_a; + int i_pdg_b; + int i_Mh; + int i_z; + int i_MX; + int i_xF; + int i_phiH; + int i_phiR; + int i_theta; + + // config options + std::set o_hadron_a_pdgs; + std::set o_hadron_b_pdgs; + std::string o_phi_r_method; + std::string o_theta_method; + enum {e_RT_via_covariant_kT} m_phi_r_method; + enum {e_hadron_a} m_theta_method; + + // storage for a single hadron + struct Hadron { + int row; + int pdg; + ROOT::Math::PxPyPzMVector p; + double z; + std::optional p_perp; + }; + + }; + +} diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml new file mode 100644 index 00000000..2469ee91 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml @@ -0,0 +1,7 @@ +physics::DihadronKinematics: + + hadron_a_list: [ 211 ] + hadron_b_list: [ -211, 211, 2212 ] + + phi_r_method: RT_via_covariant_kT + theta_method: hadron_a diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc new file mode 100644 index 00000000..a073a723 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc @@ -0,0 +1,123 @@ +#include "Validator.h" +#include "TypeDefs.h" + +#include +#include + +namespace iguana::physics { + + REGISTER_IGUANA_VALIDATOR(DihadronKinematicsValidator); + + void DihadronKinematicsValidator::Start(hipo::banklist& banks) + { + // define the algorithm sequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("physics::InclusiveKinematics"); + m_algo_seq->Add("physics::DihadronKinematics"); + m_algo_seq->SetOption("physics::DihadronKinematics", "log", m_log->GetLevel()); + m_algo_seq->SetOption>("physics::DihadronKinematics", "hadron_a_list", {particle::pi_plus}); + m_algo_seq->SetOption>("physics::DihadronKinematics", "hadron_b_list", {particle::pi_minus}); + m_algo_seq->Start(banks); + + // get bank indices + b_result = GetBankIndex(banks, "physics::DihadronKinematics"); + + // set an output file + auto output_dir = GetOutputDirectory(); + if(output_dir) { + m_output_file_basename = output_dir.value() + "/dihadron_kinematics"; + m_output_file = new TFile(m_output_file_basename + ".root", "RECREATE"); + } + + // define plots + gStyle->SetOptStat(0); + const int n_bins = 100; + plot_list = { + { + new TH1D("Mh_dist", "invariant mass M_{h} [GeV]", n_bins, 0, 4), + [](auto const& b, auto const r) { return b.getDouble("Mh", r); } + }, + { + new TH1D("z_dist", "z", n_bins, 0, 1), + [](auto const& b, auto const r) { return b.getDouble("z", r); } + }, + { + new TH1D("MX_dist", "missing mass M_X [GeV];", n_bins, 0, 4), + [](auto const& b, auto const r) { return b.getDouble("MX", r); } + }, + { + new TH1D("xF_dist", "x_{F};", n_bins, -1, 1), + [](auto const& b, auto const r) { return b.getDouble("xF", r); } + }, + { + new TH1D("phiH_dist", "#phi_{h};", n_bins, -M_PI, M_PI), + [](auto const& b, auto const r) { return b.getDouble("phiH", r); } + }, + { + new TH1D("phiR_dist", "#phi_{R}", n_bins, -M_PI, M_PI), + [](auto const& b, auto const r) { return b.getDouble("phiR", r); } + }, + { + new TH1D("theta_dist", "#theta;", n_bins, 0, M_PI), + [](auto const& b, auto const r) { return b.getDouble("theta", r); } + } + }; + + // format plots + for(auto& plot : plot_list) { + plot.hist->SetLineColor(kRed); + plot.hist->SetFillColor(kRed); + plot.hist->SetTitle( + TString(particle::title.at(particle::pi_plus)) + + TString(particle::title.at(particle::pi_minus)) + + " " + plot.hist->GetTitle()); + } + } + + + void DihadronKinematicsValidator::Run(hipo::banklist& banks) const + { + // calculate kinematics + m_algo_seq->Run(banks); + auto& result_bank = GetBank(banks, b_result, "physics::DihadronKinematics"); + + // skip events with no dihadrons + if(result_bank.getRowList().size() == 0) { + m_log->Debug("skip this event, since it has no kinematics results"); + return; + } + + // lock mutex and fill the plots + std::scoped_lock lock(m_mutex); + for(auto const& row : result_bank.getRowList()) { + for(auto& plot : plot_list) + plot.hist->Fill(plot.get_val(result_bank, row)); + } + } + + + void DihadronKinematicsValidator::Stop() + { + if(GetOutputDirectory()) { + int const n_cols = 4; + int const n_rows = (plot_list.size() - 1) / n_cols + 1; + auto canv = new TCanvas("canv", "canv", n_cols * 800, n_rows * 600); + canv->Divide(n_cols, n_rows); + int pad_num = 0; + for(auto& plot : plot_list) { + auto pad = canv->GetPad(++pad_num); + pad->cd(); + pad->SetGrid(1, 1); + pad->SetLeftMargin(0.12); + pad->SetRightMargin(0.12); + pad->SetBottomMargin(0.12); + plot.hist->Draw(); + } + canv->SaveAs(m_output_file_basename + ".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/physics/DihadronKinematics/Validator.h b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h new file mode 100644 index 00000000..d7ca5ee6 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h @@ -0,0 +1,38 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" + +#include +#include +#include +#include + +namespace iguana::physics { + + /// @brief `iguana::physics::DihadronKinematics` validator + class DihadronKinematicsValidator : public Validator + { + + DEFINE_IGUANA_VALIDATOR(DihadronKinematicsValidator, physics::DihadronKinematicsValidator) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + private: + + hipo::banklist::size_type b_result; + + struct Plot1D { + TH1D* hist; + std::function get_val; + }; + std::vector plot_list; + + TString m_output_file_basename; + TFile* m_output_file; + }; + +} diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml b/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml index 81aeabb5..144d24d1 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml @@ -46,3 +46,7 @@ actions: type: double - name: nu type: double + - name: beamPz + type: double + - name: targetM + type: double diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc index 3b52b119..5851a857 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc @@ -18,19 +18,21 @@ namespace iguana::physics { banks, b_result, GetClassName(), - {"pindex/S", "Q2/D", "x/D", "y/D", "W/D", "nu/D", "qx/D", "qy/D", "qz/D", "qE/D"}, + {"pindex/S", "Q2/D", "x/D", "y/D", "W/D", "nu/D", "qx/D", "qy/D", "qz/D", "qE/D", "beamPz/D", "targetM/D"}, 0xF000, 1); - i_pindex = result_schema.getEntryOrder("pindex"); - i_Q2 = result_schema.getEntryOrder("Q2"); - i_x = result_schema.getEntryOrder("x"); - i_y = result_schema.getEntryOrder("y"); - i_W = result_schema.getEntryOrder("W"); - i_nu = result_schema.getEntryOrder("nu"); - i_qx = result_schema.getEntryOrder("qx"); - i_qy = result_schema.getEntryOrder("qy"); - i_qz = result_schema.getEntryOrder("qz"); - i_qE = result_schema.getEntryOrder("qE"); + i_pindex = result_schema.getEntryOrder("pindex"); + i_Q2 = result_schema.getEntryOrder("Q2"); + i_x = result_schema.getEntryOrder("x"); + i_y = result_schema.getEntryOrder("y"); + i_W = result_schema.getEntryOrder("W"); + i_nu = result_schema.getEntryOrder("nu"); + i_qx = result_schema.getEntryOrder("qx"); + i_qy = result_schema.getEntryOrder("qy"); + i_qz = result_schema.getEntryOrder("qz"); + i_qE = result_schema.getEntryOrder("qE"); + i_beamPz = result_schema.getEntryOrder("beamPz"); + i_targetM = result_schema.getEntryOrder("targetM"); // instantiate RCDB reader m_rcdb = std::make_unique("RCDB|" + GetName()); @@ -111,6 +113,8 @@ namespace iguana::physics { result_bank.putDouble(i_qy, 0, result_vars.qy); result_bank.putDouble(i_qz, 0, result_vars.qz); result_bank.putDouble(i_qE, 0, result_vars.qE); + result_bank.putDouble(i_beamPz, 0, result_vars.beamPz); + result_bank.putDouble(i_targetM, 0, result_vars.targetM); ShowBank(result_bank, Logger::Header("CREATED BANK")); } @@ -244,16 +248,18 @@ namespace iguana::physics { ROOT::Math::PxPyPzMVector vec_target(target[px], target[py], target[pz], target[m]); ROOT::Math::PxPyPzMVector vec_lepton(lepton_px, lepton_py, lepton_pz, beam[m]); - auto vec_q = vec_beam - vec_lepton; - result.qx = vec_q.Px(); - result.qy = vec_q.Py(); - result.qz = vec_q.Pz(); - result.qE = vec_q.E(); - result.Q2 = -1 * vec_q.M2(); - result.x = result.Q2 / (2 * vec_q.Dot(vec_target)); - result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam); - result.W = (vec_beam + vec_target - vec_lepton).M(); - result.nu = vec_target.Dot(vec_q) / target[m]; + auto vec_q = vec_beam - vec_lepton; + result.qx = vec_q.Px(); + result.qy = vec_q.Py(); + result.qz = vec_q.Pz(); + result.qE = vec_q.E(); + result.Q2 = -1 * vec_q.M2(); + result.x = result.Q2 / (2 * vec_q.Dot(vec_target)); + result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam); + result.W = (vec_beam + vec_target - vec_lepton).M(); + result.nu = vec_target.Dot(vec_q) / target[m]; + result.beamPz = beam[pz]; + result.targetM = target[m]; return result; } diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h index dda8da34..d91a52ac 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h @@ -9,24 +9,28 @@ namespace iguana::physics { /// Set of inclusive kinematics variables struct InclusiveKinematicsVars { - /// @f$x@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{x}-component of virtual photon momentum @latex{q} vector_element_t qx; - /// @f$y@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{y}-component of virtual photon momentum @latex{q} vector_element_t qy; - /// @f$z@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{z}-component of virtual photon momentum @latex{q} vector_element_t qz; - /// @f$E@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{E}-component of virtual photon momentum @latex{q} vector_element_t qE; - /// @f$Q2@f$ (GeV@f$^2@f$) + /// @brief @latex{Q^2} (GeV@latex{^2}) double Q2; - /// @f$x_B@f$ + /// @brief @latex{x_B} double x; - /// @f$y@f$ + /// @brief @latex{y} double y; - /// @f$W@f$ (GeV) + /// @brief @latex{W} (GeV) double W; - /// @f$\nu@f$ + /// @brief @latex{\nu} double nu; + /// @brief beam momentum @latex{z}-component (GeV) + double beamPz; + /// @brief target mass (GeV) + double targetM; }; /// @brief_algo Calculate inclusive kinematics quantities defined in `iguana::physics::InclusiveKinematicsVars` @@ -63,9 +67,9 @@ namespace iguana::physics { concurrent_key_t PrepareEvent(int const runnum, double const beam_energy = -1) const; /// @action_function{scalar creator} compute kinematics from the scattered lepton. - /// @param lepton_px scattered lepton momentum component @f$p_x@f$ (GeV) - /// @param lepton_py scattered lepton momentum component @f$p_y@f$ (GeV) - /// @param lepton_pz scattered lepton momentum component @f$p_z@f$ (GeV) + /// @param lepton_px scattered lepton momentum component @latex{p_x} (GeV) + /// @param lepton_py scattered lepton momentum component @latex{p_y} (GeV) + /// @param lepton_pz scattered lepton momentum component @latex{p_z} (GeV) /// @param key the return value of `::PrepareEvent` /// @returns the reconstructed inclusive kinematics in a `iguana::physics::InclusiveKinematicsVars` instance InclusiveKinematicsVars ComputeFromLepton( @@ -104,6 +108,8 @@ namespace iguana::physics { int i_qy; int i_qz; int i_qE; + int i_beamPz; + int i_targetM; // config options mutable std::unique_ptr> o_runnum; diff --git a/src/iguana/tests/include/TestAlgorithm.h b/src/iguana/tests/include/TestAlgorithm.h index 1ccd7f06..be446b3d 100644 --- a/src/iguana/tests/include/TestAlgorithm.h +++ b/src/iguana/tests/include/TestAlgorithm.h @@ -62,15 +62,15 @@ inline int TestAlgorithm( return 1; } // print the banks, before and after - if(verbose) { - for(decltype(bank_names)::size_type it_bank = 0; it_bank < bank_names.size(); it_bank++) { - fmt::print("{:=^70}\n", fmt::format(" BEFORE: {} ", bank_names.at(it_bank))); - banks_before.at(it_bank).show(); - fmt::print("{:=^70}\n", fmt::format(" AFTER: {} ", bank_names.at(it_bank))); - banks_after.at(it_bank).show(); - fmt::print("\n"); - } - } + // if(verbose) { + // for(decltype(bank_names)::size_type it_bank = 0; it_bank < bank_names.size(); it_bank++) { + // fmt::print("{:=^70}\n", fmt::format(" BEFORE: {} ", bank_names.at(it_bank))); + // banks_before.at(it_bank).show(); + // fmt::print("{:=^70}\n", fmt::format(" AFTER: {} ", bank_names.at(it_bank))); + // banks_after.at(it_bank).show(); + // fmt::print("\n"); + // } + // } } // stop the algorithm