feat: single-hadron SIDIS kinematics (#295)

CODEOWNERS

@@ -9,10 +9,14 @@
 #################################################################################
 
 * @c-dilks
+src/iguana/algorithms/clas12/EventBuilderFilter/* @c-dilks
+src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar
+src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3
 src/iguana/algorithms/clas12/MomentumCorrection/* @RichCap @c-dilks
-src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12
+src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3
 src/iguana/algorithms/clas12/SectorFinder/* @rtysonCLAS12
-src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar
+src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12
+src/iguana/algorithms/physics/DihadronKinematics/* @c-dilks
+src/iguana/algorithms/physics/InclusiveKinematics/* @c-dilks
+src/iguana/algorithms/physics/SingleHadronKinematics/* @c-dilks
 src/iguana/services/YAMLReader.* @rtysonCLAS12 @c-dilks
-src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3
-src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3

meson/ubsan.supp

@@ -3,3 +3,4 @@ alignment:hipo::structure::getFloatAt
 alignment:hipo::structure::getShortAt
 alignment:hipo::structure::getDoubleAt
 alignment:hipo::structure::putDoubleAt
+alignment:hipo::structure::putIntAt

src/iguana/algorithms/meson.build

@@ -74,6 +74,15 @@ algo_dict = [
     'algorithm_needs_ROOT': true,
     'test_args': {'banks': [ 'REC::Particle', 'RUN::config' ]},
   },
+  {
+    'name': 'physics::SingleHadronKinematics',
+    'algorithm_needs_ROOT': true,
+    'has_action_yaml': false, # FIXME: needs a vector action function
+    'test_args': {
+      'banks': [ 'REC::Particle', 'RUN::config' ],
+      'prerequisites': [ 'physics::InclusiveKinematics' ],
+    },
+  },
   {
     'name': 'physics::DihadronKinematics',
     'algorithm_needs_ROOT': true,
@@ -92,8 +101,21 @@ algo_dict = [
 ]
 
 # make lists of objects to build; inclusion depends on whether ROOT is needed or not, and if we have ROOT
-algo_sources = [ 'Algorithm.cc', 'AlgorithmFactory.cc', 'AlgorithmSequence.cc' ]
-algo_headers = [ 'Algorithm.h', 'AlgorithmBoilerplate.h', 'TypeDefs.h', 'AlgorithmSequence.h' ]
+algo_sources = [
+  'Algorithm.cc',
+  'AlgorithmFactory.cc',
+  'AlgorithmSequence.cc',
+]
+algo_headers = [
+  'Algorithm.h',
+  'AlgorithmBoilerplate.h',
+  'TypeDefs.h',
+  'AlgorithmSequence.h',
+]
+if ROOT_dep.found()
+  algo_sources += [ 'physics/Tools.cc' ]
+  algo_headers += [ 'physics/Tools.h' ]
+endif
 vdor_sources = [ 'Validator.cc' ]
 vdor_headers = [ 'Validator.h' ]
 algo_configs = []

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

@@ -145,41 +145,41 @@ namespace iguana::physics {
       double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R());
 
       // calculate phiH
-      double phiH = PlaneAngle(
+      double phiH = tools::PlaneAngle(
           p_q.Vect(),
           p_beam.Vect(),
           p_q.Vect(),
-          p_Ph.Vect()).value_or(UNDEF);
+          p_Ph.Vect()).value_or(tools::UNDEF);
 
       // calculate PhiR
-      double phiR = UNDEF;
+      double phiR = tools::UNDEF;
       switch(m_phi_r_method) {
       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());
+            had->p_perp = tools::RejectVector(had->p.Vect(), p_q.Vect());
           }
           if(had_a.p_perp.has_value() && had_b.p_perp.has_value()) {
             auto RT = (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()) / (had_a.z + had_b.z);
-            phiR = PlaneAngle(
+            phiR = tools::PlaneAngle(
                 p_q.Vect(),
                 p_beam.Vect(),
                 p_q.Vect(),
-                RT).value_or(UNDEF);
+                RT).value_or(tools::UNDEF);
           }
           break;
         }
       }
 
       // calculate theta
-      double theta = UNDEF;
+      double theta = tools::UNDEF;
       switch(m_theta_method) {
       case e_hadron_a:
         {
-          theta = VectorAngle(
+          theta = tools::VectorAngle(
               boost__dih(had_a.p).Vect(),
-              p_Ph.Vect()).value_or(UNDEF);
+              p_Ph.Vect()).value_or(tools::UNDEF);
           break;
         }
       }
@@ -238,59 +238,6 @@ namespace iguana::physics {
 
   ///////////////////////////////////////////////////////////////////////////////
 
-  std::optional<double> DihadronKinematics::PlaneAngle(
-      ROOT::Math::XYZVector const v_a,
-      ROOT::Math::XYZVector const v_b,
-      ROOT::Math::XYZVector const v_c,
-      ROOT::Math::XYZVector const v_d)
-  {
-    auto cross_ab = v_a.Cross(v_b); // A x B
-    auto cross_cd = v_c.Cross(v_d); // C x D
-
-    auto sgn = cross_ab.Dot(v_d); // (A x B) . D
-    if(!(std::abs(sgn) > 0))
-      return std::nullopt;
-    sgn /= std::abs(sgn); // sign of (A x B) . D
-
-    auto numer = cross_ab.Dot(cross_cd); // (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 std::nullopt;
-    return sgn * std::acos(numer / denom);
-  }
-
-  std::optional<ROOT::Math::XYZVector> DihadronKinematics::ProjectVector(
-      ROOT::Math::XYZVector const v_a,
-      ROOT::Math::XYZVector const v_b)
-  {
-    auto denom = v_b.Dot(v_b);
-    if(!(std::abs(denom) > 0))
-      return std::nullopt;
-    return v_b * ( v_a.Dot(v_b) / denom );
-  }
-
-  std::optional<ROOT::Math::XYZVector> DihadronKinematics::RejectVector(
-      ROOT::Math::XYZVector const v_a,
-      ROOT::Math::XYZVector const v_b)
-  {
-    auto v_c = ProjectVector(v_a, v_b);
-    if(v_c.has_value())
-      return v_a - v_c.value();
-    return std::nullopt;
-  }
-
-  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 std::nullopt;
-  }
-
-  ///////////////////////////////////////////////////////////////////////////////
-
   void DihadronKinematics::Stop()
   {
   }

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

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "iguana/algorithms/Algorithm.h"
+#include "iguana/algorithms/physics/Tools.h"
 #include <Math/Vector3D.h>
 #include <Math/Vector4D.h>
 
@@ -18,27 +19,27 @@ namespace iguana::physics {
     int pdg_b;
     /// @brief @latex{M_h}: Invariant mass of the dihadron
     double Mh;
-    /// @brief @latex{z}: Fraction of energy of fragmenting parton carried by the dihadron
+    /// @brief @latex{z}: Momentum fraction of the fragmenting parton carried by the dihadron
     double z;
     /// @brief @latex{M_X(ehhX)}: Missing mass of the dihadron
     double MX;
     /// @brief @latex{x_F}: Feynman-x of the dihadron
     double xF;
     /// @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 `DihadronKinematics::UNDEF`, the calculation failed
+    /// if the value is `tools::UNDEF`, the calculation failed
     double phiH;
     /// @brief @latex{\phi_R}: @latex{q}-azimuthal angle between the lepton-scattering plane and dihadron plane;
-    /// if the value is `DihadronKinematics::UNDEF`, the calculation failed
+    /// if the value is `tools::UNDEF`, the calculation failed
     double phiR;
-    /// @brief @latex{\theta}: the "decay" angle of hadron A in the dihadron rest frame, with respect;
+    /// @brief @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`
   ///
   /// @begin_doc_algo{physics::DihadronKinematics | Creator}
-  /// @input_banks{REC::Particle, physics::InclusiveKinematics}
+  /// @input_banks{REC::Particle, %physics::InclusiveKinematics}
   /// @output_banks{%physics::DihadronKinematics}
   /// @end_doc
   ///
@@ -75,52 +76,11 @@ namespace iguana::physics {
       void Run(hipo::banklist& banks) const override;
       void Stop() override;
 
-      /// a value used when some calculation fails
-      double const UNDEF{-10000};
-
       /// @brief form dihadrons by pairing hadrons
       /// @param particle_bank the particle bank (`REC::Particle`)
       /// @returns a list of pairs of hadron rows
       std::vector<std::pair<int,int>> PairHadrons(hipo::bank const& particle_bank) const;
 
-      /// @brief calculate the angle between two planes
-      ///
-      /// The two planes are transverse to @latex{\vec{v}_a\times\vec{v}_b} and @latex{\vec{v}_c\times\vec{v}_d}
-      /// @param v_a vector @latex{\vec{v}_a}
-      /// @param v_b vector @latex{\vec{v}_b}
-      /// @param v_c vector @latex{\vec{v}_c}
-      /// @param v_d vector @latex{\vec{v}_d}
-      /// @returns the angle between the planes, in radians, if the calculation is successful
-      static std::optional<double> PlaneAngle(
-          ROOT::Math::XYZVector const v_a,
-          ROOT::Math::XYZVector const v_b,
-          ROOT::Math::XYZVector const v_c,
-          ROOT::Math::XYZVector const v_d);
-
-      /// @brief projection of one vector onto another
-      /// @param v_a vector @latex{\vec{v}_a}
-      /// @param v_b vector @latex{\vec{v}_b}
-      /// @returns the vector @latex{\vec{v}_a} projected onto vector @latex{\vec{v}_b}, if the calculation is successful
-      static std::optional<ROOT::Math::XYZVector> ProjectVector(
-          ROOT::Math::XYZVector const v_a,
-          ROOT::Math::XYZVector const v_b);
-
-      /// @brief projection of one vector onto the plane transverse to another vector
-      /// @param v_a vector @latex{\vec{v}_a}
-      /// @param v_b vector @latex{\vec{v}_b}
-      /// @returns the vector @latex{\vec{v}_a} projected onto the plane transverse to @latex{\vec{v}_b}, if the calculation is successful
-      static std::optional<ROOT::Math::XYZVector> RejectVector(
-          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

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

@@ -42,7 +42,7 @@ namespace iguana::physics {
         [](auto const& b, auto const r) { return b.getDouble("z", 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) { return b.getDouble("MX", r); }
       },
       {