feat: particleSvc to distribute mass by PDG

eic · Jun 4, 2024 · d8af15f · d8af15f
1 parent 2e074ae
commit d8af15f
Show file tree

Hide file tree

Showing 23 changed files with 368 additions and 26 deletions.
diff --git a/src/algorithms/interfaces/ParticleSvc.cc b/src/algorithms/interfaces/ParticleSvc.cc
diff --git a/src/algorithms/interfaces/ParticleSvc.h b/src/algorithms/interfaces/ParticleSvc.h
@@ -0,0 +1,48 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Wouter Deconinck
+
+#pragma once
+
+#include <memory>
+#include <algorithms/service.h>
+
+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>;
+
+public:
+  ParticleSvc();
+
+  virtual void init(std::shared_ptr<ParticleMap> map = nullptr) {
+    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 = 
+    std::make_shared<ParticleMap>(
+      std::make_pair(0, ParticleData{0, 0, 0.0}));
+};
+
+} // namespace algorithms
diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc
@@ -49,7 +49,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 +63,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)
       );

diff --git a/src/algorithms/reco/HadronicFinalState.h b/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
diff --git a/src/algorithms/reco/InclusiveKinematicsDA.cc b/src/algorithms/reco/InclusiveKinematicsDA.cc
@@ -47,7 +47,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 +61,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 +87,7 @@ namespace eicrecon {
     }
 
     // Calculate kinematic variables
+    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);

diff --git a/src/algorithms/reco/InclusiveKinematicsDA.h b/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
diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.cc b/src/algorithms/reco/InclusiveKinematicsElectron.cc
@@ -95,7 +95,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 +109,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 +128,7 @@ namespace eicrecon {
     }
 
     // DIS kinematics calculations
+    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);

diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.h b/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
diff --git a/src/algorithms/reco/InclusiveKinematicsJB.cc b/src/algorithms/reco/InclusiveKinematicsJB.cc
@@ -44,7 +44,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 +58,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 +75,7 @@ namespace eicrecon {
     }
 
     // Calculate kinematic variables
+    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);

diff --git a/src/algorithms/reco/InclusiveKinematicsJB.h b/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
diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.cc b/src/algorithms/reco/InclusiveKinematicsSigma.cc
@@ -46,7 +46,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)
       );
@@ -60,7 +60,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)
       );
@@ -88,6 +88,7 @@ namespace eicrecon {
     auto sigma_tot = sigma_e + sigma_h;
 
     // Calculate kinematic variables
+    const auto m_proton = m_particleSvc.particle(2122).mass;
     const auto y_sig = sigma_h / sigma_tot;
     const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
     const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);

diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.h b/src/algorithms/reco/InclusiveKinematicsSigma.h
@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm<
@@ -32,7 +34,8 @@ class InclusiveKinematicsSigma : public InclusiveKinematicsSigmaAlgorithm {
   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
diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.cc b/src/algorithms/reco/InclusiveKinematicsTruth.cc
@@ -50,7 +50,7 @@ namespace eicrecon {
     }
     const auto ei_p = ei_coll[0].getMomentum();
     const auto ei_p_mag = edm4hep::utils::magnitude(ei_p);
-    const auto ei_mass = m_electron;
+    const auto ei_mass = m_particleSvc.particle(11).mass;
     const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass));
 
     // Get incoming hadron beam
@@ -61,7 +61,7 @@ namespace eicrecon {
     }
     const auto pi_p = pi_coll[0].getMomentum();
     const auto pi_p_mag = edm4hep::utils::magnitude(pi_p);
-    const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron;
+    const auto pi_mass = m_particleSvc.particle(pi_coll[0].getPDG()).mass;
     const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass));
 
     // Get first scattered electron
@@ -76,7 +76,7 @@ namespace eicrecon {
     }
     const auto ef_p = ef_coll[0].getMomentum();
     const auto ef_p_mag = edm4hep::utils::magnitude(ef_p);
-    const auto ef_mass = m_electron;
+    const auto ef_mass = m_particleSvc.particle(11).mass;
     const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass));
 
     // DIS kinematics calculations

diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.h b/src/algorithms/reco/InclusiveKinematicsTruth.h
@@ -9,6 +9,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicsTruthAlgorithm =
@@ -29,7 +31,8 @@ class InclusiveKinematicsTruth : public InclusiveKinematicsTruthAlgorithm {
   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
diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.cc b/src/algorithms/reco/InclusiveKinematicseSigma.cc
@@ -47,7 +47,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 +61,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)
       );
@@ -96,6 +96,7 @@ namespace eicrecon {
     const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
     const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);
 
+    const auto m_proton = m_particleSvc.particle(2212).mass;
     const auto Q2_esig = Q2_e;
     const auto x_esig = x_sig;
     const auto y_esig = Q2_esig / (4.*ei.energy()*pi.energy()*x_esig); //equivalent to (2*ei.energy() / sigma_tot)*y_sig

diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.h b/src/algorithms/reco/InclusiveKinematicseSigma.h
@@ -11,6 +11,8 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
+
 namespace eicrecon {
 
 using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm<
@@ -32,7 +34,8 @@ class InclusiveKinematicseSigma : public InclusiveKinematicseSigmaAlgorithm {
   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
diff --git a/src/algorithms/reco/ScatteredElectronsTruth.cc b/src/algorithms/reco/ScatteredElectronsTruth.cc
@@ -91,12 +91,14 @@ namespace eicrecon {
       if (p.getObjectID() == ef_rc_id) {
 
         output_electrons->push_back( p );
+        static const auto m_electron = m_particleSvc.particle(11).mass;
         vScatteredElectron.SetCoordinates( p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_electron );
         electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
         // break; NOTE: if we are not computing E-Pz
         // we can safely break here and save precious CPUs
       } else {
         // Compute the sum hadronic final state
+        static const auto m_pion = m_particleSvc.particle(211).mass;
         vHadron.SetCoordinates(  p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_pion );
         vHadronicFinalState += vHadron;
       }

diff --git a/src/algorithms/reco/ScatteredElectronsTruth.h b/src/algorithms/reco/ScatteredElectronsTruth.h
@@ -10,6 +10,7 @@
 #include <string>
 #include <string_view>
 
+#include "algorithms/interfaces/ParticleSvc.h"
 
 namespace eicrecon {
 
@@ -38,7 +39,8 @@ namespace eicrecon {
     void process(const Input&, const Output&) const final;
 
   private:
-    double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_pion{0.13957};
+    const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+
   };
 
 } // namespace eicrecon