Skip to content

Commit

Permalink
feat!: SIDIS kinematics missing mass squared and Breit frame rapidity (…
Browse files Browse the repository at this point in the history 
…#308)
  • Loading branch information
@c-dilks
c-dilks authored Dec 12, 2024
1 parent 25e852b commit 81c11e2
Show file tree
Hide file tree
Showing 8 changed files with 100 additions and 31 deletions.
26 changes: 18 additions & 8 deletions src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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");
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand All @@ -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(),
Expand Down Expand Up @@ -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);
Expand Down
9 changes: 6 additions & 3 deletions src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
10 changes: 7 additions & 3 deletions src/iguana/algorithms/physics/DihadronKinematics/Validator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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); }
Expand Down
33 changes: 22 additions & 11 deletions src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
},
Expand All @@ -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");

Expand Down Expand Up @@ -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();
Expand All @@ -109,7 +114,8 @@ namespace iguana::physics {
particle_bank.getFloat("py", row),
particle_bank.getFloat("pz", row),
particle::mass.at(static_cast<particle::PDG>(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);
Expand All @@ -118,22 +124,25 @@ 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(),
p_beam.Vect(),
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);

Expand All @@ -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);
}
Expand All @@ -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);
}
Expand Down
9 changes: 6 additions & 3 deletions src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down
11 changes: 8 additions & 3 deletions src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "Validator.h"
#include "TypeDefs.h"
#include "iguana/algorithms/physics/Tools.h"

#include <Math/Vector3D.h>
#include <TStyle.h>
Expand Down Expand Up @@ -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); }
Expand Down
20 changes: 20 additions & 0 deletions src/iguana/algorithms/physics/Tools.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,5 +53,25 @@ namespace iguana::physics::tools {
return std::nullopt;
}

template <typename MOMENTUM_TYPE, typename AXIS_TYPE>
std::optional<double> 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<double> ParticleRapidity(
ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> const& momentum_vec,
ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>> const& axis_vec);
template std::optional<double> ParticleRapidity(
ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>> const& momentum_vec,
ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>> const& axis_vec);

}

13 changes: 13 additions & 0 deletions src/iguana/algorithms/physics/Tools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include <optional>
#include <Math/Vector3D.h>
#include <Math/Vector4D.h>

namespace iguana::physics::tools {

Expand Down Expand Up @@ -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 <typename MOMENTUM_TYPE, typename AXIS_TYPE>
std::optional<double> ParticleRapidity(
MOMENTUM_TYPE const& momentum_vec,
AXIS_TYPE const& axis_vec);
}

0 comments on commit 81c11e2

Please sign in to comment.