From 81c11e2803ea9e02a0b3f0912ecf64c5f9d9959a Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Thu, 12 Dec 2024 09:45:59 -0500 Subject: [PATCH] feat!: SIDIS kinematics missing mass squared and Breit frame rapidity (#308) --- .../physics/DihadronKinematics/Algorithm.cc | 26 ++++++++++----- .../physics/DihadronKinematics/Algorithm.h | 9 +++-- .../physics/DihadronKinematics/Validator.cc | 10 ++++-- .../SingleHadronKinematics/Algorithm.cc | 33 ++++++++++++------- .../SingleHadronKinematics/Algorithm.h | 9 +++-- .../SingleHadronKinematics/Validator.cc | 11 +++++-- src/iguana/algorithms/physics/Tools.cc | 20 +++++++++++ src/iguana/algorithms/physics/Tools.h | 13 ++++++++ 8 files changed, 100 insertions(+), 31 deletions(-) diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc index ecbecf05..0835c366 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc @@ -26,8 +26,9 @@ namespace iguana::physics { "Mh/D", "z/D", "PhPerp/D", - "MX/D", + "MX2/D", "xF/D", + "yB/D", "phiH/D", "phiR/D", "theta/D" @@ -41,8 +42,9 @@ namespace iguana::physics { i_Mh = result_schema.getEntryOrder("Mh"); i_z = result_schema.getEntryOrder("z"); i_PhPerp = result_schema.getEntryOrder("PhPerp"); - i_MX = result_schema.getEntryOrder("MX"); + i_MX2 = result_schema.getEntryOrder("MX2"); i_xF = result_schema.getEntryOrder("xF"); + i_yB = result_schema.getEntryOrder("yB"); i_phiH = result_schema.getEntryOrder("phiH"); i_phiR = result_schema.getEntryOrder("phiR"); i_theta = result_schema.getEntryOrder("theta"); @@ -103,11 +105,14 @@ namespace iguana::physics { inc_kin_bank.getDouble("qE", 0)); // get additional inclusive variables + auto x = inc_kin_bank.getDouble("x", 0); 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); + ROOT::Math::Boost boost__breit((p_q + 2 * x * p_target).BoostToCM()); // Breit frame + auto p_q__qp = boost__qp(p_q); + auto p_q__breit = boost__breit(p_q); // build list of dihadron rows (pindices) auto dih_rows = PairHadrons(particle_bank); @@ -130,8 +135,9 @@ namespace iguana::physics { } // calculate dihadron momenta and boosts - auto p_Ph = had_a.p + had_b.p; - auto p_Ph__qp = boost__qp(p_Ph); + auto p_Ph = had_a.p + had_b.p; + auto p_Ph__qp = boost__qp(p_Ph); + auto p_Ph__breit = boost__breit(p_Ph); ROOT::Math::Boost boost__dih(p_Ph.BoostToCM()); // CoM frame of dihadron // calculate z @@ -144,12 +150,15 @@ namespace iguana::physics { // calculate Mh double Mh = p_Ph.M(); - // calculate MX - double MX = (p_target + p_q - p_Ph).M(); + // calculate MX2 + double MX2 = (p_target + p_q - p_Ph).M2(); // calculate xF double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); + // calculate yB + double yB = tools::ParticleRapidity(p_Ph__breit, p_q__breit.Vect()).value_or(tools::UNDEF); + // calculate phiH double phiH = tools::PlaneAngle( p_q.Vect(), @@ -197,8 +206,9 @@ namespace iguana::physics { result_bank.putDouble(i_Mh, dih_row, Mh); result_bank.putDouble(i_z, dih_row, z); result_bank.putDouble(i_PhPerp, dih_row, PhPerp); - result_bank.putDouble(i_MX, dih_row, MX); + result_bank.putDouble(i_MX2, dih_row, MX2); result_bank.putDouble(i_xF, dih_row, xF); + result_bank.putDouble(i_yB, dih_row, yB); result_bank.putDouble(i_phiH, dih_row, phiH); result_bank.putDouble(i_phiR, dih_row, phiR); result_bank.putDouble(i_theta, dih_row, theta); diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h index f1d75814..40cab186 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h @@ -23,10 +23,12 @@ namespace iguana::physics { double z; /// @brief @latex{P_h^\perp}: transverse momentum of the dihadron in the @latex{\perp}-frame (transverse to @latex{\vec{q}}) double PhPerp; - /// @brief @latex{M_X(ehhX)}: Missing mass of the dihadron - double MX; + /// @brief @latex{M_X(ehhX)^2}: Missing mass squared of the dihadron + double MX2; /// @brief @latex{x_F}: Feynman-x of the dihadron double xF; + /// @brief @latex{y_{h,B}}: Breit frame rapidity of the dihadron + double yB; /// @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; @@ -98,8 +100,9 @@ namespace iguana::physics { int i_Mh; int i_z; int i_PhPerp; - int i_MX; + int i_MX2; int i_xF; + int i_yB; int i_phiH; int i_phiR; int i_theta; diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc index eeb82ab9..a8b7cdec 100644 --- a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc @@ -46,13 +46,17 @@ namespace iguana::physics { [](auto const& b, auto const r) { return b.getDouble("PhPerp", 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("MX_dist", "Missing mass: M_{X} [GeV];", n_bins, 0, 4), + [](auto const& b, auto const r) { auto MX2 = b.getDouble("MX2", r); return MX2 >= 0 ? std::sqrt(MX2) : -100; } // FIXME: handle space-like case better }, { - new TH1D("xF_dist", "x_{F};", n_bins, -1, 1), + new TH1D("xF_dist", "Feynman-x: x_{F};", n_bins, -1, 1), [](auto const& b, auto const r) { return b.getDouble("xF", r); } }, + { + new TH1D("yB_dist", "Breit frame rapidity: y_{B};", n_bins, -4, 4), + [](auto const& b, auto const r) { return b.getDouble("yB", r); } + }, { new TH1D("phiH_dist", "#phi_{h};", n_bins, -M_PI, M_PI), [](auto const& b, auto const r) { return b.getDouble("phiH", r); } diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc index 4b05acfe..47fb6fbb 100644 --- a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc @@ -24,8 +24,9 @@ namespace iguana::physics { "pdg/I", "z/D", "PhPerp/D", - "MX/D", + "MX2/D", "xF/D", + "yB/D", "phiH/D", "xi/D" }, @@ -35,8 +36,9 @@ namespace iguana::physics { i_pdg = result_schema.getEntryOrder("pdg"); i_z = result_schema.getEntryOrder("z"); i_PhPerp = result_schema.getEntryOrder("PhPerp"); - i_MX = result_schema.getEntryOrder("MX"); + i_MX2 = result_schema.getEntryOrder("MX2"); i_xF = result_schema.getEntryOrder("xF"); + i_yB = result_schema.getEntryOrder("yB"); i_phiH = result_schema.getEntryOrder("phiH"); i_xi = result_schema.getEntryOrder("xi"); @@ -81,11 +83,14 @@ namespace iguana::physics { inc_kin_bank.getDouble("qE", 0)); // get additional inclusive variables + auto x = inc_kin_bank.getDouble("x", 0); 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); + ROOT::Math::Boost boost__breit((p_q + 2 * x * p_target).BoostToCM()); // Breit frame + auto p_q__qp = boost__qp(p_q); + auto p_q__breit = boost__breit(p_q); // banks' row lists auto const& particle_bank_rowlist = particle_bank.getRowList(); @@ -109,7 +114,8 @@ namespace iguana::physics { particle_bank.getFloat("py", row), particle_bank.getFloat("pz", row), particle::mass.at(static_cast(pdg))); - auto p_Ph__qp = boost__qp(p_Ph); + auto p_Ph__qp = boost__qp(p_Ph); + auto p_Ph__breit = boost__breit(p_Ph); // calculate z double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); @@ -118,15 +124,15 @@ namespace iguana::physics { auto opt_PhPerp = tools::RejectVector(p_Ph.Vect(), p_q.Vect()); double PhPerp = opt_PhPerp.has_value() ? opt_PhPerp.value().R() : tools::UNDEF; - // 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 MX2 + double MX2 = (p_target + p_q - p_Ph).M2(); // calculate xF double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); + // calculate yB + double yB = tools::ParticleRapidity(p_Ph__breit, p_q__breit.Vect()).value_or(tools::UNDEF); + // calculate phiH double phiH = tools::PlaneAngle( p_q.Vect(), @@ -134,6 +140,9 @@ namespace iguana::physics { p_q.Vect(), p_Ph.Vect()).value_or(tools::UNDEF); + // calculate xi + double xi = p_q.Dot(p_Ph) / p_target.Dot(p_q); + // put this particle in `result_bank`'s row list result_bank_rowlist.push_back(row); @@ -142,8 +151,9 @@ namespace iguana::physics { result_bank.putInt(i_pdg, row, pdg); result_bank.putDouble(i_z, row, z); result_bank.putDouble(i_PhPerp, row, PhPerp); - result_bank.putDouble(i_MX, row, MX); + result_bank.putDouble(i_MX2, row, MX2); result_bank.putDouble(i_xF, row, xF); + result_bank.putDouble(i_yB, row, yB); result_bank.putDouble(i_phiH, row, phiH); result_bank.putDouble(i_xi, row, xi); } @@ -153,8 +163,9 @@ namespace iguana::physics { result_bank.putInt(i_pdg, row, pdg); result_bank.putDouble(i_z, row, 0); result_bank.putDouble(i_PhPerp, row, 0); - result_bank.putDouble(i_MX, row, 0); + result_bank.putDouble(i_MX2, row, 0); result_bank.putDouble(i_xF, row, 0); + result_bank.putDouble(i_yB, row, 0); result_bank.putDouble(i_phiH, row, 0); result_bank.putDouble(i_xi, row, 0); } diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h index d9fac73b..9405f5eb 100644 --- a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h @@ -16,10 +16,12 @@ namespace iguana::physics { double z; /// @brief @latex{P_h^\perp}: transverse momentum of the hadron in the @latex{\perp}-frame (transverse to @latex{\vec{q}}) double PhPerp; - /// @brief @latex{M_X(ehX)}: Missing mass of the hadron - double MX; + /// @brief @latex{M_X(ehX)^2}: Missing mass squared of the hadron + double MX2; /// @brief @latex{x_F}: Feynman-x of the hadron double xF; + /// @brief @latex{y_{h,B}}: Breit frame rapidity of the hadron + double yB; /// @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; @@ -69,8 +71,9 @@ namespace iguana::physics { int i_pdg; int i_z; int i_PhPerp; - int i_MX; + int i_MX2; int i_xF; + int i_yB; int i_phiH; int i_xi; diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc index 7c1d586e..b0c9470f 100644 --- a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc @@ -1,5 +1,6 @@ #include "Validator.h" #include "TypeDefs.h" +#include "iguana/algorithms/physics/Tools.h" #include #include @@ -41,13 +42,17 @@ namespace iguana::physics { [](auto const& b, auto const r) { return b.getDouble("PhPerp", 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("MX_dist", "Missing mass: M_{X} [GeV];", n_bins, 0, 4), + [](auto const& b, auto const r) { auto MX2 = b.getDouble("MX2", r); return MX2 >= 0 ? std::sqrt(MX2) : tools::UNDEF; } }, { - new TH1D("xF_dist", "x_{F};", n_bins, -1, 1), + new TH1D("xF_dist", "Feynman-x: x_{F};", n_bins, -1, 1), [](auto const& b, auto const r) { return b.getDouble("xF", r); } }, + { + new TH1D("yB_dist", "Breit frame rapidity: y_{B};", n_bins, -4, 4), + [](auto const& b, auto const r) { return b.getDouble("yB", r); } + }, { new TH1D("phiH_dist", "#phi_{h};", n_bins, -M_PI, M_PI), [](auto const& b, auto const r) { return b.getDouble("phiH", r); } diff --git a/src/iguana/algorithms/physics/Tools.cc b/src/iguana/algorithms/physics/Tools.cc index e06c9748..5d6166cd 100644 --- a/src/iguana/algorithms/physics/Tools.cc +++ b/src/iguana/algorithms/physics/Tools.cc @@ -53,5 +53,25 @@ namespace iguana::physics::tools { return std::nullopt; } + template + std::optional ParticleRapidity( + MOMENTUM_TYPE const& momentum_vec, + AXIS_TYPE const& axis_vec) + { + auto norm = axis_vec.R(); + if(std::abs(norm) > 0) { + auto pz = momentum_vec.Vect().Dot(axis_vec) / norm; + auto e = momentum_vec.E(); + return 0.5 * std::log((e + pz) / (e - pz)); + } + return std::nullopt; + } + template std::optional ParticleRapidity( + ROOT::Math::LorentzVector> const& momentum_vec, + ROOT::Math::DisplacementVector3D> const& axis_vec); + template std::optional ParticleRapidity( + ROOT::Math::LorentzVector> const& momentum_vec, + ROOT::Math::DisplacementVector3D> const& axis_vec); + } diff --git a/src/iguana/algorithms/physics/Tools.h b/src/iguana/algorithms/physics/Tools.h index a25686f8..77eb06d1 100644 --- a/src/iguana/algorithms/physics/Tools.h +++ b/src/iguana/algorithms/physics/Tools.h @@ -2,6 +2,7 @@ #include #include +#include namespace iguana::physics::tools { @@ -46,4 +47,16 @@ namespace iguana::physics::tools { ROOT::Math::XYZVector const v_a, ROOT::Math::XYZVector const v_b); + /// @brief calculate the rapidity of a particle, relative to an axis + /// + /// Given a particle momentum, this method calculates the rapidity + /// of the boost along an axis which takes an observer to + /// the frame in which the particle is moving perpendicular to the axis + /// @param momentum_vec the particle 4-momentum + /// @param axis_vec the axis 3-vector + /// @returns the rapidity + template + std::optional ParticleRapidity( + MOMENTUM_TYPE const& momentum_vec, + AXIS_TYPE const& axis_vec); }