diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc index 5fcf5ba1..18ce6101 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc @@ -1,6 +1,8 @@ #include "Algorithm.h" #include "TypeDefs.h" +#include + namespace iguana::physics { REGISTER_IGUANA_ALGORITHM(DihadronKinematics, "physics::DihadronKinematics"); @@ -26,7 +28,8 @@ namespace iguana::physics { "MX/D", "xF/D", "phiH/D", - "phiR/D" + "phiR/D", + "theta/D" }, 0xF000, 5); @@ -40,19 +43,27 @@ namespace iguana::physics { 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 configuration + // check phiR method if(o_phi_r_method == "RT_via_covariant_kT") - m_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)); + } /////////////////////////////////////////////////////////////////////////////// @@ -89,6 +100,13 @@ namespace iguana::physics { 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); @@ -109,13 +127,22 @@ namespace iguana::physics { particle::mass.at(static_cast(had->pdg))); } - // calculate dihadron momenta - auto p_Ph = had_a.p + had_b.p; + // 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 kinematics - double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); - double Mh = p_Ph.M(); - double MX = (p_target + p_q - p_Ph).M(); + // 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( @@ -127,22 +154,33 @@ namespace iguana::physics { // calculate PhiR double phiR = UNDEF; switch(m_phi_r_method) { - case 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 = 1 / (had_a.z + had_b.z) * (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()); - phiR = PlaneAngle( - p_q.Vect(), - p_beam.Vect(), - p_q.Vect(), - RT).value_or(UNDEF); - } - break; + 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 = 1 / (had_a.z + had_b.z) * + (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()); + 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); @@ -152,9 +190,10 @@ namespace iguana::physics { 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, 0); + 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++; } @@ -213,14 +252,12 @@ namespace iguana::physics { sgn /= std::abs(sgn); auto numer = cross_ab.Dot(cross_cd); // (A x B) . (C x D) - auto denom = std::sqrt(cross_ab.Mag2()) * std::sqrt(cross_cd.Mag2()); // |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 {}; return sgn * std::acos(numer / denom); } - /////////////////////////////////////////////////////////////////////////////// - std::optional DihadronKinematics::RejectVector( ROOT::Math::XYZVector const v_a, ROOT::Math::XYZVector const v_b) @@ -232,6 +269,16 @@ namespace iguana::physics { return v_a - v_c; } + 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 {}; + } + /////////////////////////////////////////////////////////////////////////////// void DihadronKinematics::Stop() diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h index 1baf296d..2fcfd3cb 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h @@ -30,6 +30,9 @@ namespace iguana::physics { /// @latex{\phi_R}: @latex{q}-azimuthal angle between the lepton-scattering plane and dihadron plane /// if the value is `UNDEF`, the calculation failed double phiR; + /// @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` @@ -42,7 +45,8 @@ namespace iguana::physics { /// @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 below)} + /// @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, @@ -57,6 +61,9 @@ namespace iguana::physics { /// /// @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 { @@ -95,6 +102,14 @@ namespace iguana::physics { 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 @@ -113,15 +128,15 @@ namespace iguana::physics { 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; - - enum { - RT_via_covariant_kT - } m_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 { diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml index 6e00f48a..2469ee91 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml +++ b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml @@ -4,3 +4,4 @@ physics::DihadronKinematics: 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 index c836ca41..a073a723 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc @@ -12,6 +12,7 @@ namespace iguana::physics { { // 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}); @@ -29,16 +30,47 @@ namespace iguana::physics { } // define plots - // clang-format off gStyle->SetOptStat(0); const int n_bins = 100; - Mh_dist = new TH1D("Mh_dist", "#pi^{+}#pi^{-} invariant mass M_{h};M_{h} [GeV]", n_bins, 0, 4); - // clang-format on + 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 hist : {Mh_dist}) { - hist->SetLineColor(kRed); - hist->SetFillColor(kRed); + 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()); } } @@ -49,6 +81,7 @@ namespace iguana::physics { 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; @@ -57,7 +90,8 @@ namespace iguana::physics { // lock mutex and fill the plots std::scoped_lock lock(m_mutex); for(auto const& row : result_bank.getRowList()) { - Mh_dist->Fill(result_bank.getDouble("Mh", row)); + for(auto& plot : plot_list) + plot.hist->Fill(plot.get_val(result_bank, row)); } } @@ -65,22 +99,19 @@ namespace iguana::physics { void DihadronKinematicsValidator::Stop() { if(GetOutputDirectory()) { - int n_rows = 2; - int n_cols = 4; + 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); - for(int pad_num = 1; pad_num <= n_rows * n_cols; pad_num++) { - auto pad = canv->GetPad(pad_num); + 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); - switch(pad_num) { - case 1: - Mh_dist->Draw(); - break; - } + plot.hist->Draw(); } canv->SaveAs(m_output_file_basename + ".png"); m_output_file->Write(); diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.h b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h index 12442564..d7ca5ee6 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Validator.h +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h @@ -25,7 +25,11 @@ namespace iguana::physics { hipo::banklist::size_type b_result; - mutable TH1D* Mh_dist; + struct Plot1D { + TH1D* hist; + std::function get_val; + }; + std::vector plot_list; TString m_output_file_basename; TFile* m_output_file;