From 79c50ea6e54919fe64a8caec8b0b1304466b7024 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Wed, 23 Oct 2024 19:22:27 -0400 Subject: [PATCH] feat: single-hadron SIDIS kinematics (#295) --- CODEOWNERS | 12 +- meson/ubsan.supp | 1 + src/iguana/algorithms/meson.build | 26 ++- .../physics/DihadronKinematics/Algorithm.cc | 71 +------- .../physics/DihadronKinematics/Algorithm.h | 52 +----- .../physics/DihadronKinematics/Validator.cc | 2 +- .../SingleHadronKinematics/Algorithm.cc | 167 ++++++++++++++++++ .../SingleHadronKinematics/Algorithm.h | 79 +++++++++ .../SingleHadronKinematics/Config.yaml | 3 + .../SingleHadronKinematics/Validator.cc | 113 ++++++++++++ .../SingleHadronKinematics/Validator.h | 38 ++++ src/iguana/algorithms/physics/Tools.cc | 57 ++++++ src/iguana/algorithms/physics/Tools.h | 49 +++++ 13 files changed, 555 insertions(+), 115 deletions(-) create mode 100644 src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc create mode 100644 src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h create mode 100644 src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml create mode 100644 src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc create mode 100644 src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h create mode 100644 src/iguana/algorithms/physics/Tools.cc create mode 100644 src/iguana/algorithms/physics/Tools.h diff --git a/CODEOWNERS b/CODEOWNERS index 565ee6e9..1c777213 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -9,10 +9,14 @@ ################################################################################# * @c-dilks +src/iguana/algorithms/clas12/EventBuilderFilter/* @c-dilks +src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar +src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3 src/iguana/algorithms/clas12/MomentumCorrection/* @RichCap @c-dilks -src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12 +src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3 src/iguana/algorithms/clas12/SectorFinder/* @rtysonCLAS12 -src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar +src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12 +src/iguana/algorithms/physics/DihadronKinematics/* @c-dilks +src/iguana/algorithms/physics/InclusiveKinematics/* @c-dilks +src/iguana/algorithms/physics/SingleHadronKinematics/* @c-dilks src/iguana/services/YAMLReader.* @rtysonCLAS12 @c-dilks -src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3 -src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3 \ No newline at end of file diff --git a/meson/ubsan.supp b/meson/ubsan.supp index fe8381cd..51719faa 100644 --- a/meson/ubsan.supp +++ b/meson/ubsan.supp @@ -3,3 +3,4 @@ alignment:hipo::structure::getFloatAt alignment:hipo::structure::getShortAt alignment:hipo::structure::getDoubleAt alignment:hipo::structure::putDoubleAt +alignment:hipo::structure::putIntAt diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 49331342..7ac67bcb 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::SingleHadronKinematics', + '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': 'physics::DihadronKinematics', 'algorithm_needs_ROOT': true, @@ -92,8 +101,21 @@ algo_dict = [ ] # make lists of objects to build; inclusion depends on whether ROOT is needed or not, and if we have ROOT -algo_sources = [ 'Algorithm.cc', 'AlgorithmFactory.cc', 'AlgorithmSequence.cc' ] -algo_headers = [ 'Algorithm.h', 'AlgorithmBoilerplate.h', 'TypeDefs.h', 'AlgorithmSequence.h' ] +algo_sources = [ + 'Algorithm.cc', + 'AlgorithmFactory.cc', + 'AlgorithmSequence.cc', +] +algo_headers = [ + 'Algorithm.h', + 'AlgorithmBoilerplate.h', + 'TypeDefs.h', + 'AlgorithmSequence.h', +] +if ROOT_dep.found() + algo_sources += [ 'physics/Tools.cc' ] + algo_headers += [ 'physics/Tools.h' ] +endif vdor_sources = [ 'Validator.cc' ] vdor_headers = [ 'Validator.h' ] algo_configs = [] diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc index 6b8a95b4..0705d988 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc @@ -145,41 +145,41 @@ namespace iguana::physics { double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); // calculate phiH - double phiH = PlaneAngle( + double phiH = tools::PlaneAngle( p_q.Vect(), p_beam.Vect(), p_q.Vect(), - p_Ph.Vect()).value_or(UNDEF); + p_Ph.Vect()).value_or(tools::UNDEF); // calculate PhiR - double phiR = UNDEF; + double phiR = tools::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()); + had->p_perp = tools::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( + phiR = tools::PlaneAngle( p_q.Vect(), p_beam.Vect(), p_q.Vect(), - RT).value_or(UNDEF); + RT).value_or(tools::UNDEF); } break; } } // calculate theta - double theta = UNDEF; + double theta = tools::UNDEF; switch(m_theta_method) { case e_hadron_a: { - theta = VectorAngle( + theta = tools::VectorAngle( boost__dih(had_a.p).Vect(), - p_Ph.Vect()).value_or(UNDEF); + p_Ph.Vect()).value_or(tools::UNDEF); break; } } @@ -238,59 +238,6 @@ namespace iguana::physics { /////////////////////////////////////////////////////////////////////////////// - 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 index 2a761b2e..ca8dd510 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h @@ -1,6 +1,7 @@ #pragma once #include "iguana/algorithms/Algorithm.h" +#include "iguana/algorithms/physics/Tools.h" #include #include @@ -18,19 +19,19 @@ namespace iguana::physics { 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 + /// @brief @latex{z}: Momentum fraction of the 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 + /// if the value is `tools::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 + /// if the value is `tools::UNDEF`, the calculation failed double phiR; - /// @brief @latex{\theta}: the "decay" angle of hadron A in the dihadron rest frame, with respect; + /// @brief @latex{\theta}: The "decay" angle of hadron A in the dihadron rest frame, with respect; /// to the dihadron momentum direction double theta; }; @@ -38,7 +39,7 @@ namespace iguana::physics { /// @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} + /// @input_banks{REC::Particle, %physics::InclusiveKinematics} /// @output_banks{%physics::DihadronKinematics} /// @end_doc /// @@ -75,52 +76,11 @@ namespace iguana::physics { 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 diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc index a073a723..6e936c6d 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc @@ -42,7 +42,7 @@ namespace iguana::physics { [](auto const& b, auto const r) { return b.getDouble("z", r); } }, { - new TH1D("MX_dist", "missing mass M_X [GeV];", n_bins, 0, 4), + new TH1D("MX_dist", "missing mass M_{X} [GeV];", n_bins, 0, 4), [](auto const& b, auto const r) { return b.getDouble("MX", r); } }, { diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc new file mode 100644 index 00000000..d0fb6f4c --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc @@ -0,0 +1,167 @@ +#include "Algorithm.h" +#include "TypeDefs.h" +#include "iguana/algorithms/physics/Tools.h" + +#include + +namespace iguana::physics { + + REGISTER_IGUANA_ALGORITHM(SingleHadronKinematics, "physics::SingleHadronKinematics"); + + void SingleHadronKinematics::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/S", + "pdg/I", + "z/D", + "MX/D", + "xF/D", + "phiH/D", + "xi/D" + }, + 0xF000, + 7); + i_pindex = result_schema.getEntryOrder("pindex"); + i_pdg = result_schema.getEntryOrder("pdg"); + 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_xi = result_schema.getEntryOrder("xi"); + + // parse config file + ParseYAMLConfig(); + o_hadron_pdgs = GetOptionSet("hadron_list"); + + } + + /////////////////////////////////////////////////////////////////////////////// + + void SingleHadronKinematics::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); + + // banks' row lists + auto const& particle_bank_rowlist = particle_bank.getRowList(); + hipo::bank::rowlist::list_t result_bank_rowlist{}; + result_bank.setRows(particle_bank.getRows()); + + // loop over ALL rows of `particle_bank` + // - we will calculate kinematics for rows in `particle_bank_rowlist`, and zero out all the other rows + // - we want the `result_bank` to have the same number of rows as `particle_bank` and the same ordering, + // so that banks which reference `particle_bank` rows can be used to reference `result_bank` rows too + for(int row = 0; row < particle_bank.getRows(); row++) { + + // if the particle is in `o_hadron_pdgs` AND the row is in `particle_bank`'s filtered row list + if(auto pdg{particle_bank.getInt("pid", row)}; + o_hadron_pdgs.find(pdg) != o_hadron_pdgs.end() && + std::find(particle_bank_rowlist.begin(), particle_bank_rowlist.end(), row) != particle_bank_rowlist.end()) { + + // hadron momentum + auto p_Ph = ROOT::Math::PxPyPzMVector( + particle_bank.getFloat("px", row), + particle_bank.getFloat("py", row), + particle_bank.getFloat("pz", row), + particle::mass.at(static_cast(pdg))); + auto p_Ph__qp = boost__qp(p_Ph); + + // calculate z + double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); + + // calculate xi + double xi = p_q.Dot(p_Ph) / p_target.Dot(p_q); + + // 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 = tools::PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + p_Ph.Vect()).value_or(tools::UNDEF); + + // put this particle in `result_bank`'s row list + result_bank_rowlist.push_back(row); + + // fill the bank + result_bank.putShort(i_pindex, row, row); + result_bank.putInt(i_pdg, row, pdg); + result_bank.putDouble(i_z, row, z); + result_bank.putDouble(i_MX, row, MX); + result_bank.putDouble(i_xF, row, xF); + result_bank.putDouble(i_phiH, row, phiH); + result_bank.putDouble(i_xi, row, xi); + } + else { + // zero the row + result_bank.putShort(i_pindex, row, row); + result_bank.putInt(i_pdg, row, pdg); + result_bank.putDouble(i_z, row, 0); + result_bank.putDouble(i_MX, row, 0); + result_bank.putDouble(i_xF, row, 0); + result_bank.putDouble(i_phiH, row, 0); + result_bank.putDouble(i_xi, row, 0); + } + } + + // apply the filtered rowlist to `result_bank` + result_bank.getMutableRowList().setList(result_bank_rowlist); + + ShowBank(result_bank, Logger::Header("CREATED BANK")); + } + + /////////////////////////////////////////////////////////////////////////////// + + void SingleHadronKinematics::Stop() + { + } + +} diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h new file mode 100644 index 00000000..f98d9adb --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h @@ -0,0 +1,79 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" +#include +#include + +namespace iguana::physics { + + /// Set of hadron kinematics variables + struct SingleHadronKinematicsVars { + /// @brief `REC::Particle` row (`pindex`) of the hadron + int pindex; + /// @brief PDG code of the hadron + int pdg; + /// @brief @latex{z}: Momentum fraction of the fragmenting parton carried by the hadron + double z; + /// @brief @latex{M_X(ehX)}: Missing mass of the hadron + double MX; + /// @brief @latex{x_F}: Feynman-x of the hadron + 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 `tools::UNDEF`, the calculation failed + double phiH; + /// @brief @latex{\xi_h}: Longitudinal momentum fraction of the nucleon carried by the hadron + double xi; + }; + + /// @brief_algo Calculate semi-inclusive hadron kinematic quantities defined in `iguana::physics::SingleHadronKinematicsVars` + /// + /// @begin_doc_algo{physics::SingleHadronKinematics | Creator} + /// @input_banks{REC::Particle, %physics::InclusiveKinematics} + /// @output_banks{%physics::SingleHadronKinematics} + /// @end_doc + /// + /// @begin_doc_config + /// @config_param{hadron_list | list[int] | list of hadron PDGs} + /// @end_doc + /// + /// The output bank `%physics::SingleHadronKinematics` will have the same number of rows as the input particle bank `REC::Particle` + /// - we want the output bank to have the same number of rows and ordering as the input + /// particle bank, so that banks which reference the input particle bank's rows (usually via `pindex`) can be used to + /// reference the output bank's rows too + /// - rows of the input particle bank which were filtered upstream will also be filtered out here, and all the values of the + /// corresponding row in the output bank will be zeroed, since no calculations are performed for + /// those particles + /// - particles which are not listed in the configuration parameter `hadron_list` will also be filtered out and zeroed + class SingleHadronKinematics : public Algorithm + { + + DEFINE_IGUANA_ALGORITHM(SingleHadronKinematics, physics::SingleHadronKinematics) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + 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; + int i_pdg; + int i_z; + int i_MX; + int i_xF; + int i_phiH; + int i_xi; + + // config options + std::set o_hadron_pdgs; + + }; + +} diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml b/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml new file mode 100644 index 00000000..b2277af3 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml @@ -0,0 +1,3 @@ +physics::SingleHadronKinematics: + + hadron_list: [ 211, -211, 2212 ] diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc new file mode 100644 index 00000000..d84edeb9 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc @@ -0,0 +1,113 @@ +#include "Validator.h" +#include "TypeDefs.h" + +#include +#include + +namespace iguana::physics { + + REGISTER_IGUANA_VALIDATOR(SingleHadronKinematicsValidator); + + void SingleHadronKinematicsValidator::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::SingleHadronKinematics"); + m_algo_seq->SetOption("physics::SingleHadronKinematics", "log", m_log->GetLevel()); + m_algo_seq->SetOption>("physics::SingleHadronKinematics", "hadron_list", {particle::pi_plus}); + m_algo_seq->Start(banks); + + // get bank indices + b_result = GetBankIndex(banks, "physics::SingleHadronKinematics"); + + // set an output file + auto output_dir = GetOutputDirectory(); + if(output_dir) { + m_output_file_basename = output_dir.value() + "/single_hadron_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("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("xi_dist", "#xi", n_bins, -1, 1), + [](auto const& b, auto const r) { return b.getDouble("xi", r); } + }, + }; + + // format plots + for(auto& plot : plot_list) { + plot.hist->SetLineColor(kGreen+1); + plot.hist->SetFillColor(kGreen+1); + plot.hist->SetTitle( + TString(particle::title.at(particle::pi_plus)) + + " " + plot.hist->GetTitle()); + } + } + + + void SingleHadronKinematicsValidator::Run(hipo::banklist& banks) const + { + // calculate kinematics + m_algo_seq->Run(banks); + auto& result_bank = GetBank(banks, b_result, "physics::SingleHadronKinematics"); + + // skip events with no hadrons + 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 SingleHadronKinematicsValidator::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/SingleHadronKinematics/Validator.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h new file mode 100644 index 00000000..3017da57 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h @@ -0,0 +1,38 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" + +#include +#include +#include +#include + +namespace iguana::physics { + + /// @brief `iguana::physics::SingleHadronKinematics` validator + class SingleHadronKinematicsValidator : public Validator + { + + DEFINE_IGUANA_VALIDATOR(SingleHadronKinematicsValidator, physics::SingleHadronKinematicsValidator) + + 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/Tools.cc b/src/iguana/algorithms/physics/Tools.cc new file mode 100644 index 00000000..e06c9748 --- /dev/null +++ b/src/iguana/algorithms/physics/Tools.cc @@ -0,0 +1,57 @@ +#include "Tools.h" + +namespace iguana::physics::tools { + + 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) + { + 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 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 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 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; + } + +} + diff --git a/src/iguana/algorithms/physics/Tools.h b/src/iguana/algorithms/physics/Tools.h new file mode 100644 index 00000000..a25686f8 --- /dev/null +++ b/src/iguana/algorithms/physics/Tools.h @@ -0,0 +1,49 @@ +/// @file Tools.h + +#include +#include + +namespace iguana::physics::tools { + + /// a value used when some calculation fails + double const UNDEF{-10000}; + + /// @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 + 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 + 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 + 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 + std::optional VectorAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + +}