feat: dihadron rapidity

JeffersonLab · Nov 20, 2024 · 8a2095a · 8a2095a
1 parent 3b66440
commit 8a2095a
Show file tree

Hide file tree

Showing 5 changed files with 24 additions and 7 deletions.
diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc
@@ -28,6 +28,7 @@ namespace iguana::physics {
           "PhPerp/D",
           "MX2/D",
           "xF/D",
+          "yB/D",
           "phiH/D",
           "phiR/D",
           "theta/D"
@@ -43,6 +44,7 @@ namespace iguana::physics {
     i_PhPerp   = result_schema.getEntryOrder("PhPerp");
     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
@@ -150,6 +156,9 @@ namespace iguana::physics {
       // 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(),
@@ -199,6 +208,7 @@ namespace iguana::physics {
       result_bank.putDouble(i_PhPerp,  dih_row, PhPerp);
       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
@@ -27,6 +27,8 @@ namespace iguana::physics {
     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;
@@ -100,6 +102,7 @@ namespace iguana::physics {
       int i_PhPerp;
       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
@@ -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),
+        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.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h
@@ -20,7 +20,7 @@ namespace iguana::physics {
     double MX2;
     /// @brief @latex{x_F}: Feynman-x of the hadron
     double xF;
-    /// @brief @latex{y_{h,B}}: Breit frame rapidity
+    /// @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

diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc
@@ -42,7 +42,7 @@ 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),
+        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; }
       },
       {