Merge remote-tracking branch 'origin/main' into uniq-new-banks

.github/workflows/ci.yml

@@ -389,7 +389,7 @@ jobs:
             -Dtest_data_file=$(pwd)/test_data.hipo       \
             -Dtest_num_events=${{ env.num_events }}      \
             -Dtest_output_dir=$(pwd)/validation_results  \
-            -Dtest_multithreading=${{ env.num_threads }} \
+            -Dtest_num_threads=${{ env.num_threads }}    \
             ${{ matrix.opts }}
       - name: dump all build options
         run: meson configure iguana_build --no-pager
@@ -413,6 +413,7 @@ jobs:
       - run: tree iguana
       ### run tests
       - name: meson test
+        if: ${{ matrix.id != 'documentation' }}
         working-directory: iguana_build
         run: |
           if [ "${{ env.verbose_test }}" = "true" ]; then
@@ -423,26 +424,26 @@ jobs:
           fi
       ### run benchmarks
       - name: benchmark algorithms
-        if: ${{ matrix.id == 'coverage' }} # use the coverage job's GITHUB_STEP_SUMMARY
+        if: ${{ matrix.id == 'cpp' && inputs.id == 'linux-latest' }}
         run: |
           for suite in single_threaded memoize; do
             meson test --benchmark --suite $suite -C iguana_build | tee benchmark_$suite.txt
           done
           iguana_src/.github/make-benchmark-table.rb benchmark_{single_threaded,memoize}.txt | xargs -0 -I{} echo {} >> $GITHUB_STEP_SUMMARY
       ### coverage
-      # - name: coverage
-      #   if: ${{ matrix.id == 'coverage' }}
-      #   run: |
-      #     ninja -C iguana_build coverage-html
-      #     ninja -C iguana_build coverage-text
-      #     mv iguana_build/meson-logs/coveragereport coverage-report
-      #     echo '### Coverage Report' >> $GITHUB_STEP_SUMMARY
-      #     echo '```' >> $GITHUB_STEP_SUMMARY
-      #     cat iguana_build/meson-logs/coverage.txt >> $GITHUB_STEP_SUMMARY
-      #     echo '```' >> $GITHUB_STEP_SUMMARY
-      #     echo '' >> $GITHUB_STEP_SUMMARY
-      #     echo '- for details, see the `coverage-report` artifact' >> $GITHUB_STEP_SUMMARY
-      #     echo '- to compare to the report from the `main` branch, see <https://jeffersonlab.github.io/iguana/coverage-report>' >> $GITHUB_STEP_SUMMARY
+      - name: coverage
+        if: ${{ matrix.id == 'coverage' }}
+        run: |
+          ninja -C iguana_build coverage-html
+          ninja -C iguana_build coverage-text
+          mv iguana_build/meson-logs/coveragereport coverage-report
+          echo '### Coverage Report' >> $GITHUB_STEP_SUMMARY
+          echo '```' >> $GITHUB_STEP_SUMMARY
+          cat iguana_build/meson-logs/coverage.txt >> $GITHUB_STEP_SUMMARY
+          echo '```' >> $GITHUB_STEP_SUMMARY
+          echo '' >> $GITHUB_STEP_SUMMARY
+          echo '- for details, see the `coverage-report` artifact' >> $GITHUB_STEP_SUMMARY
+          echo '- to compare to the report from the `main` branch, see <https://jeffersonlab.github.io/iguana/coverage-report>' >> $GITHUB_STEP_SUMMARY
       ### test relocatability
       - name: test relocatability
         if: ${{ matrix.id == 'cpp' }} # don't bother re-running santizers, etc.
@@ -511,13 +512,13 @@ jobs:
           name: meson_logs_${{ matrix.id }}
           retention-days: 5
           path: iguana_build/meson-logs
-      # - name: upload coverage artifacts
-      #   uses: actions/upload-artifact@v4
-      #   if: ${{ matrix.id == 'coverage' }}
-      #   with:
-      #     name: coverage-report
-      #     retention-days: 5
-      #     path: coverage-report
+      - name: upload coverage artifacts
+        uses: actions/upload-artifact@v4
+        if: ${{ matrix.id == 'coverage' }}
+        with:
+          name: coverage-report
+          retention-days: 5
+          path: coverage-report
       - name: upload validator artifacts
         uses: actions/upload-artifact@v4
         if: ${{ matrix.id == 'cpp' }}
@@ -551,18 +552,18 @@ jobs:
         with:
           name: doc_doxygen
           path: doxygen
-      # - name: download coverage report
-      #   uses: actions/download-artifact@v4
-      #   with:
-      #     name: coverage-report
-      #     path: coverage-report
+      - name: download coverage report
+        uses: actions/download-artifact@v4
+        with:
+          name: coverage-report
+          path: coverage-report
       - run: tree
       - name: collect
         run: |
           mkdir pages
           cp iguana_src/.github/pages-index.html pages/index.html
           mv doxygen pages/
-      # - run: mv coverage-report pages/
+      - run: mv coverage-report pages/
       - run: tree
       - uses: actions/upload-pages-artifact@v3
         with:

.github/workflows/linux.yml

@@ -36,7 +36,7 @@ jobs:
           ],
           "include": [
             { "id": "cpp",                 "CC": "gcc",   "CXX": "g++",     "opts": "-Dbuildtype=release -Dz_require_root=true   -Dtest_validator_all_stats=true"  },
-            { "id": "coverage",            "CC": "gcc",   "CXX": "g++",     "opts": "-Dbuildtype=release -Dz_require_root=true   -Db_coverage=true"     },
+            { "id": "coverage",            "CC": "gcc",   "CXX": "g++",     "opts": "-Dbuildtype=release -Dz_require_root=true   -Db_coverage=true  -Dz_test_multithreading=false" },
             { "id": "documentation",       "CC": "gcc",   "CXX": "g++",     "opts": "-Dbuildtype=release -Dz_require_root=true   -Dinstall_documentation=true" },
             { "id": "address-sanitizer",   "CC": "clang", "CXX": "clang++", "opts": "-Dbuildtype=debug   -Dz_require_root=true   -Db_sanitize=address   -Db_lundef=false -Db_pie=true" },
             { "id": "thread-sanitizer",    "CC": "clang", "CXX": "clang++", "opts": "-Dbuildtype=debug   -Dz_require_root=true   -Db_sanitize=thread    -Db_lundef=false -Db_pie=true" },

meson.options

@@ -11,9 +11,10 @@ option('test_data_file',           type: 'string',  value: '',    description: '
 option('test_num_events',          type: 'string',  value: '10',  description: 'Number of events from `test_data_file` to test')
 option('test_output_dir',          type: 'string',  value: '',    description: 'Output directory for tests. Must be an absolute path. If unspecified, tests will still run, but will not produce output files.')
 option('test_validator_all_stats', type: 'boolean', value: false, description: 'If true, use all statistics for validators, rather than `test_num_events`')
-option('test_multithreading',      type: 'integer', value: 4, min: 0, description: 'run multithreading tests with this many threads (use 0 for hardware max)')
+option('test_num_threads',         type: 'integer', value: 4, min: 0, description: 'run multithreading tests with this many threads (use 0 for hardware max)')
 
 # expert options: the defaults should be reasonable for a local installation; different values may be preferred for installation in common areas
 option('z_install_envfile',     type: 'boolean', value: true,  description: 'Install a sourceable environment variable file')
 option('z_require_root',        type: 'boolean', value: false, description: 'Fail if ROOT is not found')
 option('z_require_rcdb',        type: 'boolean', value: false, description: 'Fail if RCDB is not found')
+option('z_test_multithreading', type: 'boolean', value: true,  description: 'Enable multithreading tests')

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

@@ -21,8 +21,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");
@@ -83,11 +84,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);
@@ -110,8 +114,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
@@ -124,12 +129,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(),
@@ -177,8 +185,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

@@ -67,8 +67,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

@@ -19,8 +19,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");
 
@@ -65,11 +66,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();
@@ -93,7 +97,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);
@@ -102,22 +107,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);
 
@@ -126,8 +134,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);
       }
@@ -137,8 +146,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

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