feat!: SIDIS kinematics missing mass squared and Breit frame rapidity

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);

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;

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); }

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<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);
@@ -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);
 
@@ -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);
       }

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;
 

src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc

@@ -1,5 +1,6 @@
 #include "Validator.h"
 #include "TypeDefs.h"
+#include "iguana/algorithms/physics/Tools.h"
 
 #include <Math/Vector3D.h>
 #include <TStyle.h>
@@ -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); }

src/iguana/algorithms/physics/Tools.cc

@@ -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);
+
 }
 

src/iguana/algorithms/physics/Tools.h

@@ -2,6 +2,7 @@
 
 #include <optional>
 #include <Math/Vector3D.h>
+#include <Math/Vector4D.h>
 
 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 <typename MOMENTUM_TYPE, typename AXIS_TYPE>
+  std::optional<double> ParticleRapidity(
+      MOMENTUM_TYPE const& momentum_vec,
+      AXIS_TYPE const& axis_vec);
 }