Reverted where the proton filter is selected from

eic · Aug 15, 2024 · feaa087 · feaa087
1 parent ce9e4e4
commit feaa087
Show file tree

Hide file tree

Showing 5 changed files with 25 additions and 12 deletions.
diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc
@@ -34,14 +34,13 @@ void eicrecon::MatrixTransferStatic::process(
     const MatrixTransferStatic::Input& input,
     const MatrixTransferStatic::Output& output) {
 
-  const auto [beamP,scatP,rechits] = input;
+  const auto [mcparts,rechits] = input;
   auto [outputParticles] = output;
 
   //Set beam energy from first MCBeamElectron, using std::call_once
   std::call_once(m_initBeamE,[&](){
 
-    m_initialized = initalizeMatrix(*beamP);
-    if(!m_initialized) m_initialized = initalizeMatrix(*scatP);
+    m_initialized = initalizeMatrix(*mcparts);
 
   });
 
@@ -170,12 +169,27 @@ bool eicrecon::MatrixTransferStatic::initalizeMatrix(const edm4hep::MCParticleCo
     return false;
   }
 
+  double numBeamProtons = 0;
+  double runningMomentum = 0.0;
+
+  for (const auto p: mcparts) {
+    if(mcparts.size() == 1 && p.getPDG() == 2212){
+      runningMomentum = p.getMomentum().z;
+      numBeamProtons++;
+    }
+    if (p.getGeneratorStatus() == 4 && p.getPDG() == 2212) { //look for "beam" proton
+      runningMomentum += p.getMomentum().z;
+      numBeamProtons++;
+    }
+  }
+
+  if(numBeamProtons == 0) {error("No beam protons to choose matrix!"); return false;}
+
+  double nomMomentum = runningMomentum/numBeamProtons;
+
   // Fractional error allowed in beam momentum
   double nomMomentumError = 0.05;
 
-  // Set nominal momentum to MCParticle momentum
-  double nomMomentum = mcparts.at(0).getMomentum().z;
-
   double aX[2][2] = {{0.0, 0.0},
                     {0.0, 0.0}};
   double aY[2][2] = {{0.0, 0.0},

diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.h b/src/algorithms/fardetectors/MatrixTransferStatic.h
@@ -13,6 +13,7 @@
 #include <gsl/pointers>
 #include <string>
 #include <string_view>
+#include <Eigen/Dense>
 
 #include "MatrixTransferStaticConfig.h"
 #include "algorithms/interfaces/WithPodConfig.h"
@@ -21,7 +22,6 @@ namespace eicrecon {
 
   using MatrixTransferStaticAlgorithm = algorithms::Algorithm<
     algorithms::Input<
-      edm4hep::MCParticleCollection,
       edm4hep::MCParticleCollection,
       edm4eic::TrackerHitCollection
     >,

diff --git a/src/detectors/FOFFMTRK/FOFFMTRK.cc b/src/detectors/FOFFMTRK/FOFFMTRK.cc
@@ -69,7 +69,7 @@ void InitPlugin(JApplication *app) {
 
     app->Add(new JOmniFactoryGeneratorT<MatrixTransferStatic_factory>(
         "ForwardOffMRecParticles",
-        {"MCBeamProtons","MCScatteredProtons","ForwardOffMTrackerRecHits"},
+        {"MCParticles","ForwardOffMTrackerRecHits"},
         {"ForwardOffMRecParticles"},
         recon_cfg,
         app

diff --git a/src/detectors/RPOTS/RPOTS.cc b/src/detectors/RPOTS/RPOTS.cc
@@ -70,7 +70,7 @@ void InitPlugin(JApplication *app) {
 
     app->Add(new JOmniFactoryGeneratorT<MatrixTransferStatic_factory>(
         "ForwardRomanPotRecParticles",
-        {"MCBeamProtons","MCScatteredProtons","ForwardRomanPotRecHits"},
+        {"MCParticles","ForwardRomanPotRecHits"},
         {"ForwardRomanPotRecParticles"},
         recon_cfg,
         app

diff --git a/src/factories/fardetectors/MatrixTransferStatic_factory.h b/src/factories/fardetectors/MatrixTransferStatic_factory.h
@@ -25,8 +25,7 @@ class MatrixTransferStatic_factory :
 private:
     std::unique_ptr<AlgoT> m_algo;
 
-    PodioInput<edm4hep::MCParticle> m_mcbeamp_input {this};
-    PodioInput<edm4hep::MCParticle> m_mcscatp_input {this};
+    PodioInput<edm4hep::MCParticle> m_mcparts_input {this};
     PodioInput<edm4eic::TrackerHit> m_hits_input {this};
     PodioOutput<edm4eic::ReconstructedParticle> m_tracks_output {this};
 
@@ -68,7 +67,7 @@ class MatrixTransferStatic_factory :
     }
 
     void Process(int64_t run_number, uint64_t event_number) {
-        m_algo->process({m_mcbeamp_input(), m_mcscatp_input(), m_hits_input()}, {m_tracks_output().get()});
+        m_algo->process({m_mcparts_input(), m_hits_input()}, {m_tracks_output().get()});
     }
 
 };