feat: add more variables, and plot all of them

src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc

@@ -1,6 +1,8 @@
 #include "Algorithm.h"
 #include "TypeDefs.h"
 
+#include <Math/Boost.h>
+
 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<int>("hadron_a_list");
     o_hadron_b_pdgs = GetOptionSet<int>("hadron_b_list");
     o_phi_r_method  = GetOptionScalar<std::string>("phi_r_method");
+    o_theta_method  = GetOptionScalar<std::string>("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<particle::PDG>(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<ROOT::Math::XYZVector> 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<double> 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()

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<double> 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<int> o_hadron_a_pdgs;
       std::set<int> 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 {

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

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

@@ -12,6 +12,7 @@ namespace iguana::physics {
   {
     // define the algorithm sequence
     m_algo_seq = std::make_unique<AlgorithmSequence>();
+    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<std::vector<int>>("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,30 +90,28 @@ namespace iguana::physics {
     // lock mutex and fill the plots
     std::scoped_lock<std::mutex> 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));
     }
   }
 
 
   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();

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<double(hipo::bank const&, int const)> get_val;
+      };
+      std::vector<Plot1D> plot_list;
 
       TString m_output_file_basename;
       TFile* m_output_file;