feat: particleSvc to distribute mass by PDG (#1487)

### Briefly, what does this PR introduce?
Instead of storing the mass of the proton etc in various classes, this
PR adds an algorithms::ParticleSvc that distributes the correct masses
and charge etc by PDG number.

Note: I was playing around with using the algorithms::ParticleSvc as a
thin interface-only, and let a JService read a data file and provide the
actual service. But I gave up since it seemed pointless.

### What kind of change does this PR introduce?
- [ ] Bug fix (issue #__)
- [x] New feature (issue #__)
- [ ] 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?
No.

### Does this PR change default behavior?
No.

---------

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

src/algorithms/interfaces/CMakeLists.txt

@@ -1,2 +1,10 @@
 set(PLUGIN_NAME "algorithms_interfaces")
-plugin_headers_only(${PLUGIN_NAME})
+
+# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets
+# Setting default includes, libraries and installation paths
+plugin_add(${PLUGIN_NAME} WITH_SHARED_LIBRARY WITHOUT_PLUGIN)
+
+# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then
+# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds
+# headers to the correct installation directory
+plugin_glob_all(${PLUGIN_NAME})

src/algorithms/interfaces/ParticleSvc.h

@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Wouter Deconinck
+
+#pragma once
+
+#include <algorithms/service.h>
+#include <map>
+#include <memory>
+#include <string>
+
+namespace algorithms {
+
+class ParticleSvc : public Service<ParticleSvc> {
+public:
+  struct ParticleData {
+    int         pdgCode;
+    int         charge;
+    double      mass;
+    std::string name;
+  };
+  using Particle    = ParticleData;
+  using ParticleMap = std::map<int, Particle>;
+
+private:
+  static const std::shared_ptr<ParticleMap> kParticleMap;
+
+public:
+  virtual void init(std::shared_ptr<ParticleMap> map = kParticleMap) {
+    if (map != nullptr) {
+      m_particleMap = map;
+    }
+  }
+
+  virtual std::shared_ptr<ParticleMap> particleMap() const {
+    return m_particleMap;
+  };
+
+  virtual Particle& particle(int pdg) const {
+    if (m_particleMap->count(pdg) == 0) {
+      return m_particleMap->at(0);
+    }
+    return m_particleMap->at(pdg);
+  };
+
+protected:
+  std::shared_ptr<ParticleMap> m_particleMap{nullptr};
+
+  ALGORITHMS_DEFINE_SERVICE(ParticleSvc)
+};
+
+} // namespace algorithms

src/algorithms/pid/IrtCherenkovParticleID.cc

@@ -111,11 +111,9 @@ void IrtCherenkovParticleID::init(
   }
 
   // get PDG info for the particles we want to identify in PID
-  // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we
-  // have a proper PDG database service, we hard-code the masses in Tools.h
   m_log->debug("List of particles for PID:");
   for(auto pdg : m_cfg.pdgList) {
-    auto mass = Tools::GetPDGMass(pdg);
+    auto mass = m_particleSvc.particle(pdg).mass;
     m_pdg_mass.insert({ pdg, mass });
     m_log->debug("  {:>8}  M={} GeV", pdg, mass);
   }

src/algorithms/pid/IrtCherenkovParticleID.h

@@ -21,6 +21,7 @@
 
 // EICrecon
 #include "IrtCherenkovParticleIDConfig.h"
+#include "algorithms/interfaces/ParticleSvc.h"
 #include "algorithms/interfaces/WithPodConfig.h"
 
 namespace eicrecon {
@@ -69,6 +70,8 @@ namespace eicrecon {
     CherenkovDetectorCollection*    m_irt_det_coll;
     CherenkovDetector*              m_irt_det;
 
+    const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+
     uint64_t    m_cell_mask;
     std::string m_det_name;
     std::unordered_map<int,double>           m_pdg_mass;

src/algorithms/pid/Tools.h

@@ -62,34 +62,6 @@ namespace eicrecon {
       }
 
 
-      // -------------------------------------------------------------------------------------
-      // PDG mass lookup
-
-      // local PDG mass database
-      // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we
-      // have a proper PDG database service, we hard-code the masses we need;
-      // use Tools::GetPDGMass for access
-      static std::unordered_map<int,double> GetPDGMasses() {
-        return std::unordered_map<int,double>{
-          { 11,   0.000510999 },
-          { 211,  0.13957     },
-          { 321,  0.493677    },
-          { 2212, 0.938272    }
-        };
-      }
-
-      static double GetPDGMass(int pdg) {
-        double mass;
-        try { mass = GetPDGMasses().at(std::abs(pdg)); }
-        catch(const std::out_of_range& e) {
-          throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown PDG={} in algorithms/pid/Tools::GetPDGMass",pdg));
-        }
-        return mass;
-      }
-
-      static int GetNumPDGs() { return GetPDGMasses().size(); };
-
-
       // -------------------------------------------------------------------------------------
       // Table rebinning and lookup
 

src/algorithms/reco/HadronicFinalState.cc

@@ -25,13 +25,7 @@ using ROOT::Math::PxPyPzEVector;
 
 namespace eicrecon {
 
-  void HadronicFinalState::init() {
-    // m_pidSvc = service("ParticleSvc");
-    // if (!m_pidSvc) {
-    //   debug("Unable to locate Particle Service. "
-    //     "Make sure you have ParticleSvc in the configuration.");
-    // }
-  }
+  void HadronicFinalState::init() { }
 
   void HadronicFinalState::process(
       const HadronicFinalState::Input& input,
@@ -49,7 +43,7 @@ namespace eicrecon {
     const PxPyPzEVector ei(
       round_beam_four_momentum(
         ei_coll[0].getMomentum(),
-        m_electron,
+        m_particleSvc.particle(ei_coll[0].getPDG()).mass,
         {-5.0, -10.0, -18.0},
         0.0)
       );
@@ -63,7 +57,7 @@ namespace eicrecon {
     const PxPyPzEVector pi(
       round_beam_four_momentum(
         pi_coll[0].getMomentum(),
-        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
+        m_particleSvc.particle(pi_coll[0].getPDG()).mass,
         {41.0, 100.0, 275.0},
         m_crossingAngle)
       );

src/algorithms/reco/HadronicFinalState.h

@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using HadronicFinalStateAlgorithm = algorithms::Algorithm<
@@ -31,7 +33,8 @@ class HadronicFinalState : public HadronicFinalStateAlgorithm {
   void process(const Input&, const Output&) const final;
 
 private:
-  double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025};
+  const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+  double m_crossingAngle{-0.025};
 };
 
 } // namespace eicrecon

src/algorithms/reco/InclusiveKinematicsDA.cc

@@ -23,13 +23,7 @@ using ROOT::Math::PxPyPzEVector;
 
 namespace eicrecon {
 
-  void InclusiveKinematicsDA::init() {
-    // m_pidSvc = service("ParticleSvc");
-    // if (!m_pidSvc) {
-    //   m_log->debug("Unable to locate Particle Service. "
-    //     "Make sure you have ParticleSvc in the configuration.");
-    // }
-  }
+  void InclusiveKinematicsDA::init() { }
 
   void InclusiveKinematicsDA::process(
       const InclusiveKinematicsDA::Input& input,
@@ -47,7 +41,7 @@ namespace eicrecon {
     const PxPyPzEVector ei(
       round_beam_four_momentum(
         ei_coll[0].getMomentum(),
-        m_electron,
+        m_particleSvc.particle(ei_coll[0].getPDG()).mass,
         {-5.0, -10.0, -18.0},
         0.0)
       );
@@ -61,7 +55,7 @@ namespace eicrecon {
     const PxPyPzEVector pi(
       round_beam_four_momentum(
         pi_coll[0].getMomentum(),
-        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
+        m_particleSvc.particle(pi_coll[0].getPDG()).mass,
         {41.0, 100.0, 275.0},
         m_crossingAngle)
       );
@@ -87,6 +81,7 @@ namespace eicrecon {
     }
 
     // Calculate kinematic variables
+    static const auto m_proton = m_particleSvc.particle(2212).mass;
     const auto y_da = tan(gamma_h/2.) / ( tan(theta_e/2.) + tan(gamma_h/2.) );
     const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(gamma_h/2.)) );
     const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);

src/algorithms/reco/InclusiveKinematicsDA.h

@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm<
@@ -32,7 +34,8 @@ class InclusiveKinematicsDA : public InclusiveKinematicsDAAlgorithm {
   void process(const Input&, const Output&) const final;
 
 private:
-  double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025};
+  const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+  double m_crossingAngle{-0.025};
 };
 
 } // namespace eicrecon

src/algorithms/reco/InclusiveKinematicsElectron.cc

@@ -23,13 +23,7 @@ using ROOT::Math::PxPyPzEVector;
 
 namespace eicrecon {
 
-  void InclusiveKinematicsElectron::init() {
-    // m_pidSvc = service("ParticleSvc");
-    // if (!m_pidSvc) {
-    //   debug("Unable to locate Particle Service. "
-    //     "Make sure you have ParticleSvc in the configuration.");
-    // }
-  }
+  void InclusiveKinematicsElectron::init() { }
 
   void InclusiveKinematicsElectron::process(
       const InclusiveKinematicsElectron::Input& input,
@@ -95,7 +89,7 @@ namespace eicrecon {
     const PxPyPzEVector ei(
       round_beam_four_momentum(
         ei_coll[0].getMomentum(),
-        m_electron,
+        m_particleSvc.particle(ei_coll[0].getPDG()).mass,
         {-5.0, -10.0, -18.0},
         0.0)
       );
@@ -109,7 +103,7 @@ namespace eicrecon {
     const PxPyPzEVector pi(
       round_beam_four_momentum(
         pi_coll[0].getMomentum(),
-        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
+        m_particleSvc.particle(pi_coll[0].getPDG()).mass,
         {41.0, 100.0, 275.0},
         m_crossingAngle)
       );
@@ -128,6 +122,7 @@ namespace eicrecon {
     }
 
     // DIS kinematics calculations
+    static const auto m_proton = m_particleSvc.particle(2212).mass;
     const auto ef = electrons.front();
     const auto q = ei - ef;
     const auto q_dot_pi = q.Dot(pi);

src/algorithms/reco/InclusiveKinematicsElectron.h

@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm<
@@ -32,7 +34,8 @@ class InclusiveKinematicsElectron : public InclusiveKinematicsElectronAlgorithm
   void process(const Input&, const Output&) const final;
 
 private:
-  double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025};
+  const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+  double m_crossingAngle{-0.025};
 };
 
 } // namespace eicrecon

src/algorithms/reco/InclusiveKinematicsJB.cc

@@ -20,13 +20,7 @@ using ROOT::Math::PxPyPzEVector;
 
 namespace eicrecon {
 
-  void InclusiveKinematicsJB::init() {
-    // m_pidSvc = service("ParticleSvc");
-    // if (!m_pidSvc) {
-    //   debug("Unable to locate Particle Service. "
-    //     "Make sure you have ParticleSvc in the configuration.");
-    // }
-  }
+  void InclusiveKinematicsJB::init() { }
 
   void InclusiveKinematicsJB::process(
       const InclusiveKinematicsJB::Input& input,
@@ -44,7 +38,7 @@ namespace eicrecon {
     const PxPyPzEVector ei(
       round_beam_four_momentum(
         ei_coll[0].getMomentum(),
-        m_electron,
+        m_particleSvc.particle(ei_coll[0].getPDG()).mass,
         {-5.0, -10.0, -18.0},
         0.0)
       );
@@ -58,7 +52,7 @@ namespace eicrecon {
     const PxPyPzEVector pi(
       round_beam_four_momentum(
         pi_coll[0].getMomentum(),
-        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
+        m_particleSvc.particle(pi_coll[0].getPDG()).mass,
         {41.0, 100.0, 275.0},
         m_crossingAngle)
       );
@@ -75,6 +69,7 @@ namespace eicrecon {
     }
 
     // Calculate kinematic variables
+    static const auto m_proton = m_particleSvc.particle(2212).mass;
     const auto y_jb = sigma_h / (2.*ei.energy());
     const auto Q2_jb = ptsum*ptsum / (1. - y_jb);
     const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb);

src/algorithms/reco/InclusiveKinematicsJB.h

@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm<
@@ -32,7 +34,8 @@ class InclusiveKinematicsJB : public InclusiveKinematicsJBAlgorithm {
   void process(const Input&, const Output&) const final;
 
 private:
-  double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025};
+  const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+  double m_crossingAngle{-0.025};
 };
 
 } // namespace eicrecon