Created algorithm and factory for calculating hadronic final state va…

…riables (#1405)

### Briefly, what does this PR introduce?

This adds and algorithm and corresponding factory for calculating
hadronic final state variables. Since this uses a new `EDM4eic` data
type not included in major release 5, the algorithm is only called for
major release >= 6.

### What kind of change does this PR introduce?
- [ ] Bug fix (issue #__)
- [x] New feature (issue #1362 )
- [ ] Documentation update
- [ ] Other: __

### Please check if this PR fulfills the following:
- [ ] Tests for the changes have been added
- [ ] Documentation has been added / updated
- [ ] Changes have been communicated to collaborators

### Does this PR introduce breaking changes? What changes might users
need to make to their code?

### Does this PR change default behavior?

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Dmitry Kalinkin <[email protected]>

src/algorithms/reco/HadronicFinalState.cc

@@ -0,0 +1,154 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Tyler Kutz
+
+#include <edm4eic/EDM4eicVersion.h>
+#if EDM4EIC_VERSION_MAJOR >= 6
+
+#include <Math/GenVector/LorentzVector.h>
+#include <Math/GenVector/PxPyPzE4D.h>
+#include <Math/Vector4Dfwd.h>
+#include <edm4eic/HadronicFinalStateCollection.h>
+#include <edm4eic/MCRecoParticleAssociationCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/Vector3f.h>
+#include <fmt/core.h>
+#include <podio/ObjectID.h>
+#include <cmath>
+#include <gsl/pointers>
+#include <vector>
+
+#include "Beam.h"
+#include "Boost.h"
+#include "HadronicFinalState.h"
+
+using ROOT::Math::PxPyPzEVector;
+
+namespace eicrecon {
+
+  void HadronicFinalState::init(std::shared_ptr<spdlog::logger>& logger) {
+    m_log = logger;
+    // m_pidSvc = service("ParticleSvc");
+    // if (!m_pidSvc) {
+    //   m_log->debug("Unable to locate Particle Service. "
+    //     "Make sure you have ParticleSvc in the configuration.");
+    // }
+  }
+
+  void HadronicFinalState::process(
+      const HadronicFinalState::Input& input,
+      const HadronicFinalState::Output& output) const {
+
+    const auto [mcparts, rcparts, rcassoc] = input;
+    auto [hadronicfinalstate] = output;
+
+    // Get incoming electron beam
+    const auto ei_coll = find_first_beam_electron(mcparts);
+    if (ei_coll.size() == 0) {
+      m_log->debug("No beam electron found");
+      return;
+    }
+    const PxPyPzEVector ei(
+      round_beam_four_momentum(
+        ei_coll[0].getMomentum(),
+        m_electron,
+        {-5.0, -10.0, -18.0},
+        0.0)
+      );
+
+    // Get incoming hadron beam
+    const auto pi_coll = find_first_beam_hadron(mcparts);
+    if (pi_coll.size() == 0) {
+      m_log->debug("No beam hadron found");
+      return;
+    }
+    const PxPyPzEVector pi(
+      round_beam_four_momentum(
+        pi_coll[0].getMomentum(),
+        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
+        {41.0, 100.0, 275.0},
+        m_crossingAngle)
+      );
+
+    // Get first scattered electron
+    const auto ef_coll = find_first_scattered_electron(mcparts);
+    if (ef_coll.size() == 0) {
+      m_log->debug("No truth scattered electron found");
+      return;
+    }
+    // Associate first scattered electron with reconstructed electrons
+    //const auto ef_assoc = std::find_if(
+    //  rcassoc->begin(),
+    //  rcassoc->end(),
+    //  [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); });
+    auto ef_assoc = rcassoc->begin();
+    for (; ef_assoc != rcassoc->end(); ++ef_assoc) {
+      if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) {
+        break;
+      }
+    }
+    if (!(ef_assoc != rcassoc->end())) {
+      m_log->debug("Truth scattered electron not in reconstructed particles");
+      return;
+    }
+    const auto ef_rc{ef_assoc->getRec()};
+    const auto ef_rc_id{ef_rc.getObjectID().index};
+
+    // Sums in colinear frame
+    double pxsum = 0;
+    double pysum = 0;
+    double pzsum = 0;
+    double Esum = 0;
+
+    // Get boost to colinear frame
+    auto boost = determine_boost(ei, pi);
+
+    auto hfs = hadronicfinalstate->create(0., 0., 0.);
+
+    for (const auto& p: *rcparts) {
+
+      bool isHadron = true;
+      // Check if it's the scattered electron
+      if (p.getObjectID().index == ef_rc_id) isHadron = false;
+      // Check for non-hadron PDG codes
+      if (p.getPDG() == 11) isHadron = false;
+      if (p.getPDG() == 22) isHadron = false;
+      if (p.getPDG() == 13) isHadron = false;
+      // If it's the scattered electron or not a hadron, skip
+      if(!isHadron) continue;
+
+      // Lorentz vector in lab frame
+      PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
+      // Boost to colinear frame
+      PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab);
+
+      pxsum += hf_boosted.Px();
+      pysum += hf_boosted.Py();
+      pzsum += hf_boosted.Pz();
+      Esum += hf_boosted.E();
+
+      hfs.addToHadrons(p);
+
+    }
+
+    // Hadronic final state calculations
+    auto sigma = Esum - pzsum;
+    auto pT = sqrt(pxsum*pxsum + pysum*pysum);
+    auto gamma = (pT*pT - sigma*sigma)/(pT*pT + sigma*sigma);
+
+    hfs.setSigma(sigma);
+    hfs.setPT(pT);
+    hfs.setGamma(gamma);
+
+    // Sigma zero or negative
+    if (sigma <= 0) {
+      m_log->debug("Sigma zero or negative");
+      return;
+    }
+
+    m_log->debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma());
+
+  }
+
+} // namespace Jug::Reco
+#endif

src/algorithms/reco/HadronicFinalState.h

@@ -0,0 +1,48 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Tyler Kutz
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/HadronicFinalStateCollection.h>
+#include <edm4eic/MCRecoParticleAssociationCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <spdlog/logger.h>
+#include <memory>
+#include <string>
+#include <string_view>
+
+
+namespace eicrecon {
+
+  using HadronicFinalStateAlgorithm = algorithms::Algorithm<
+    algorithms::Input<
+      edm4hep::MCParticleCollection,
+      edm4eic::ReconstructedParticleCollection,
+      edm4eic::MCRecoParticleAssociationCollection
+    >,
+    algorithms::Output<
+      edm4eic::HadronicFinalStateCollection
+    >
+  >;
+
+  class HadronicFinalState
+  : public HadronicFinalStateAlgorithm {
+
+  public:
+    HadronicFinalState(std::string_view name)
+      : HadronicFinalStateAlgorithm{name,
+                            {"MCParticles", "inputParticles", "inputAssociations"},
+                            {"hadronicFinalState"},
+                            "Calculate summed quantities of the hadronic final state."} {}
+
+    void init(std::shared_ptr<spdlog::logger>& logger);
+    void process(const Input&, const Output&) const final;
+
+  private:
+    std::shared_ptr<spdlog::logger> m_log;
+    double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025};
+  };
+
+} // namespace eicrecon

src/factories/reco/HadronicFinalState_factory.h

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Tyler Kutz
+
+#pragma once
+
+#include <JANA/JEvent.h>
+#include <edm4eic/HadronicFinalStateCollection.h>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "extensions/jana/JOmniFactory.h"
+
+namespace eicrecon {
+
+template <typename AlgoT>
+class HadronicFinalState_factory :
+        public JOmniFactory<HadronicFinalState_factory<AlgoT>> {
+
+public:
+    using FactoryT = JOmniFactory<HadronicFinalState_factory<AlgoT>>;
+
+private:
+    std::unique_ptr<AlgoT> m_algo;
+
+    typename FactoryT::template PodioInput<edm4hep::MCParticle> m_mc_particles_input {this};
+    typename FactoryT::template PodioInput<edm4eic::ReconstructedParticle> m_rc_particles_input {this};
+    typename FactoryT::template PodioInput<edm4eic::MCRecoParticleAssociation> m_rc_particles_assoc_input {this};
+    typename FactoryT::template PodioOutput<edm4eic::HadronicFinalState> m_hadronic_final_state_output {this};
+
+public:
+    void Configure() {
+        m_algo = std::make_unique<AlgoT>(this->GetPrefix());
+        m_algo->init(this->logger());
+    }
+
+    void ChangeRun(int64_t run_number) {
+    }
+
+    void Process(int64_t run_number, uint64_t event_number) {
+        m_algo->process({m_mc_particles_input(), m_rc_particles_input(), m_rc_particles_assoc_input()},
+                        {m_hadronic_final_state_output().get()});
+    }
+};
+
+} // eicrecon

src/global/reco/reco.cc

@@ -8,6 +8,7 @@
 #include <edm4eic/MCRecoClusterParticleAssociation.h>
 #include <edm4eic/MCRecoParticleAssociation.h>
 #include <edm4eic/ReconstructedParticle.h>
+#include <edm4eic/EDM4eicVersion.h>
 #include <edm4hep/MCParticle.h>
 #include <algorithm>
 #include <map>
@@ -19,6 +20,9 @@
 #include "algorithms/reco/InclusiveKinematicsJB.h"
 #include "algorithms/reco/InclusiveKinematicsSigma.h"
 #include "algorithms/reco/InclusiveKinematicseSigma.h"
+#if EDM4EIC_VERSION_MAJOR >= 6
+#include "algorithms/reco/HadronicFinalState.h"
+#endif
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
 #include "factories/meta/CollectionCollector_factory.h"
 #include "factories/meta/FilterMatching_factory.h"
@@ -27,6 +31,9 @@
 #include "factories/reco/InclusiveKinematicsTruth_factory.h"
 #include "factories/reco/JetReconstruction_factory.h"
 #include "factories/reco/TransformBreitFrame_factory.h"
+#if EDM4EIC_VERSION_MAJOR >= 6
+#include "factories/reco/HadronicFinalState_factory.h"
+#endif
 #include "global/reco/ChargedReconstructedParticleSelector_factory.h"
 #include "global/reco/MC2SmearedParticle_factory.h"
 #include "global/reco/MatchClusters_factory.h"
@@ -277,5 +284,20 @@ void InitPlugin(JApplication *app) {
             app
     ));
 
+#if EDM4EIC_VERSION_MAJOR >= 6
+    app->Add(new JOmniFactoryGeneratorT<HadronicFinalState_factory<HadronicFinalState>>(
+        "HadronicFinalState",
+        {
+          "MCParticles",
+          "ReconstructedParticles",
+          "ReconstructedParticleAssociations"
+        },
+        {
+          "HadronicFinalState"
+        },
+        app
+    ));
+#endif
+
 }
 } // extern "C"

src/services/io/podio/JEventProcessorPODIO.cc

@@ -1,6 +1,8 @@
 
 #include "JEventProcessorPODIO.h"
 
+#include <edm4eic/EDM4eicVersion.h>
+
 #include <JANA/JApplication.h>
 #include <JANA/JLogger.h>
 #include <JANA/Services/JParameterManager.h>
@@ -182,6 +184,9 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             "ReconstructedElectrons",
             "ScatteredElectronsTruth",
             "ScatteredElectronsEMinusPz",
+#if EDM4EIC_VERSION_MAJOR >= 6
+            "HadronicFinalState",
+#endif
 
             // Track projections
             "CalorimeterTrackProjections",