From 04807aea1d821b6dd0adfe461048da8e85045649 Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 13 Aug 2024 13:00:27 +0100 Subject: [PATCH 1/6] Updated MatrixTransferStatic to only initalize once and use pre-filtered MCparticle collections --- .../fardetectors/MatrixTransferStatic.cc | 279 ++++++++++-------- .../fardetectors/MatrixTransferStatic.h | 16 +- src/detectors/FOFFMTRK/FOFFMTRK.cc | 8 +- src/detectors/RPOTS/RPOTS.cc | 8 +- .../MatrixTransferStatic_factory.h | 5 +- 5 files changed, 184 insertions(+), 132 deletions(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index 48331e6c67..9319072e3a 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.cc +++ b/src/algorithms/fardetectors/MatrixTransferStatic.cc @@ -26,151 +26,42 @@ void eicrecon::MatrixTransferStatic::init() { void eicrecon::MatrixTransferStatic::process( const MatrixTransferStatic::Input& input, - const MatrixTransferStatic::Output& output) const { + const MatrixTransferStatic::Output& output) { - const auto [mcparts, rechits] = input; + const auto [beamP,scatP,rechits] = input; auto [outputParticles] = output; std::vector> aX(m_cfg.aX); std::vector> aY(m_cfg.aY); //----- Define constants here ------ - double aXinv[2][2] = {{0.0, 0.0}, + double m_aXinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}}; - double aYinv[2][2] = {{0.0, 0.0}, + double m_aYinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}}; - double nomMomentum = m_cfg.nomMomentum; //extract the nominal value first -- will be overwritten by MCParticle - double local_x_offset = m_cfg.local_x_offset; - double local_y_offset = m_cfg.local_y_offset; + double nomMomentum = m_cfg.nomMomentum; //extract the nominal value first -- will be overwritten by MCParticle + double local_x_offset = m_cfg.local_x_offset; + double local_y_offset = m_cfg.local_y_offset; double local_x_slope_offset = m_cfg.local_x_slope_offset; double local_y_slope_offset = m_cfg.local_y_slope_offset; 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!! Skipping!!"); return;} - - nomMomentum = runningMomentum/numBeamProtons; - - double nomMomentumError = 0.05; - - //This is a temporary solution to get the beam energy information - //needed to select the correct matrix - - if(abs(275.0 - nomMomentum)/275.0 < nomMomentumError){ - - aX[0][0] = 3.251116; //a - aX[0][1] = 30.285734; //b - aX[1][0] = 0.186036375; //c - aX[1][1] = 0.196439472; //d - - aY[0][0] = 0.4730500000; //a - aY[0][1] = 3.062999454; //b - aY[1][0] = 0.0204108951; //c - aY[1][1] = -0.139318692; //d - - local_x_offset = -0.339334; - local_y_offset = -0.000299454; - local_x_slope_offset = -0.219603248; - local_y_slope_offset = -0.000176128; - - } - else if(abs(100.0 - nomMomentum)/100.0 < nomMomentumError){ - - aX[0][0] = 3.152158; //a - aX[0][1] = 20.852072; //b - aX[1][0] = 0.181649517; //c - aX[1][1] = -0.303998487; //d - - aY[0][0] = 0.5306100000; //a - aY[0][1] = 3.19623343; //b - aY[1][0] = 0.0226283320; //c - aY[1][1] = -0.082666019; //d - - local_x_offset = -0.329072; - local_y_offset = -0.00028343; - local_x_slope_offset = -0.218525084; - local_y_slope_offset = -0.00015321; - - } - else if(abs(41.0 - nomMomentum)/41.0 < nomMomentumError){ + //Set beam energy from first MCBeamElectron, using std::call_once + std::call_once(m_initBeamE,[&](){ - aX[0][0] = 3.135997; //a - aX[0][1] = 18.482273; //b - aX[1][0] = 0.176479921; //c - aX[1][1] = -0.497839483; //d + m_initialized = initalizeMatrix(*beamP); + if(!m_initialized) m_initialized = initalizeMatrix(*scatP); - aY[0][0] = 0.4914400000; //a - aY[0][1] = 4.53857451; //b - aY[1][0] = 0.0179664765; //c - aY[1][1] = 0.004160679; //d + }); - local_x_offset = -0.283273; - local_y_offset = -0.00552451; - local_x_slope_offset = -0.21174031; - local_y_slope_offset = -0.003212011; - - } - else if(abs(135.0 - nomMomentum)/135.0 < nomMomentumError){ //135 GeV deuterons - - aX[0][0] = 1.6248; - aX[0][1] = 12.966293; - aX[1][0] = 0.1832; - aX[1][1] = -2.8636535; - - aY[0][0] = 0.0001674; //a - aY[0][1] = -28.6003; //b - aY[1][0] = 0.0000837; //c - aY[1][1] = -2.87985; //d - - local_x_offset = -11.9872; - local_y_offset = -0.0146; - local_x_slope_offset = -14.75315; - local_y_slope_offset = -0.0073; - - } - else { - error("MatrixTransferStatic:: No valid matrix found to match beam momentum!! Skipping!!"); + if(!m_initialized){ + debug("No valid matrix can be identified for reconstruction"); return; } - - double determinant = aX[0][0] * aX[1][1] - aX[0][1] * aX[1][0]; - - if (determinant == 0) { - error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!"); - return; - } - - aXinv[0][0] = aX[1][1] / determinant; - aXinv[0][1] = -aX[0][1] / determinant; - aXinv[1][0] = -aX[1][0] / determinant; - aXinv[1][1] = aX[0][0] / determinant; - - - determinant = aY[0][0] * aY[1][1] - aY[0][1] * aY[1][0]; - - if (determinant == 0) { - error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!"); - return; - } - - aYinv[0][0] = aY[1][1] / determinant; - aYinv[0][1] = -aY[0][1] / determinant; - aYinv[1][0] = -aY[1][0] / determinant; - aYinv[1][1] = aY[0][0] / determinant; - + //---- begin Reconstruction code ---- edm4hep::Vector3f goodHit[2] = {{0.0,0.0,0.0},{0.0,0.0,0.0}}; @@ -247,8 +138,8 @@ void eicrecon::MatrixTransferStatic::process( for (unsigned i0 = 0; i0 < 2; i0++) { for (unsigned i1 = 0; i1 < 2; i1++) { - Xip[i0] += aXinv[i0][i1] * Xrp[i1]; - Yip[i0] += aYinv[i0][i1] * Yrp[i1]; + Xip[i0] += m_aXinv[i0][i1] * Xrp[i1]; + Yip[i0] += m_aYinv[i0][i1] * Yrp[i1]; } } @@ -257,7 +148,7 @@ void eicrecon::MatrixTransferStatic::process( double rsy = Yip[1] * dd4hep::mrad; // calculate momentum magnitude from measured deltaP – using thin lens optics. - double p = nomMomentum * (1 + 0.01 * Xip[0]); + double p = m_nomMomentum * (1 + 0.01 * Xip[0]); double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy); edm4hep::Vector3f prec = {static_cast(p * rsx / norm), static_cast(p * rsy / norm), @@ -281,3 +172,137 @@ void eicrecon::MatrixTransferStatic::process( } // end enough hits for matrix reco } + +// Initalize the matrix based on the MC particle collections +bool eicrecon::MatrixTransferStatic::initalizeMatrix(const edm4hep::MCParticleCollection& mcparts) { + + // Check size of MC particle collection + if(mcparts.size() == 0){ + error("No MCParticles found in collection!"); + return false; + } + + // 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}, + {0.0, 0.0}}; + + //This is a temporary solution to get the beam energy information + //needed to select the correct matrix + if(abs(275.0 - nomMomentum)/275.0 < nomMomentumError){ + + aX[0][0] = 3.251116; //a + aX[0][1] = 30.285734; //b + aX[1][0] = 0.186036375; //c + aX[1][1] = 0.196439472; //d + + aY[0][0] = 0.4730500000; //a + aY[0][1] = 3.062999454; //b + aY[1][0] = 0.0204108951; //c + aY[1][1] = -0.139318692; //d + + m_local_x_offset = -0.339334; + m_local_y_offset = -0.000299454; + m_local_x_slope_offset = -0.219603248; + m_local_y_slope_offset = -0.000176128; + + m_nomMomentum = 275.0; + + } + else if(abs(100.0 - nomMomentum)/100.0 < nomMomentumError){ + + aX[0][0] = 3.152158; //a + aX[0][1] = 20.852072; //b + aX[1][0] = 0.181649517; //c + aX[1][1] = -0.303998487; //d + + aY[0][0] = 0.5306100000; //a + aY[0][1] = 3.19623343; //b + aY[1][0] = 0.0226283320; //c + aY[1][1] = -0.082666019; //d + + m_local_x_offset = -0.329072; + m_local_y_offset = -0.00028343; + m_local_x_slope_offset = -0.218525084; + m_local_y_slope_offset = -0.00015321; + + m_nomMomentum = 100.0; + + } + else if(abs(41.0 - nomMomentum)/41.0 < nomMomentumError){ + + aX[0][0] = 3.135997; //a + aX[0][1] = 18.482273; //b + aX[1][0] = 0.176479921; //c + aX[1][1] = -0.497839483; //d + + aY[0][0] = 0.4914400000; //a + aY[0][1] = 4.53857451; //b + aY[1][0] = 0.0179664765; //c + aY[1][1] = 0.004160679; //d + + m_local_x_offset = -0.283273; + m_local_y_offset = -0.00552451; + m_local_x_slope_offset = -0.21174031; + m_local_y_slope_offset = -0.003212011; + + m_nomMomentum = 41.0; + + } + else if(abs(135.0 - nomMomentum)/135.0 < nomMomentumError){ //135 GeV deuterons + + aX[0][0] = 1.6248; + aX[0][1] = 12.966293; + aX[1][0] = 0.1832; + aX[1][1] = -2.8636535; + + aY[0][0] = 0.0001674; //a + aY[0][1] = -28.6003; //b + aY[1][0] = 0.0000837; //c + aY[1][1] = -2.87985; //d + + m_local_x_offset = -11.9872; + m_local_y_offset = -0.0146; + m_local_x_slope_offset = -14.75315; + m_local_y_slope_offset = -0.0073; + + m_nomMomentum = 135.0; + + } + else { + error("No valid matrix found to match beam momentum!"); + return false; + } + + double determinantX = aX[0][0] * aX[1][1] - aX[0][1] * aX[1][0]; + + if (determinantX == 0) { + error("Reco x matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!"); + return false; + } + + m_aXinv[0][0] = aX[1][1] / determinantX; + m_aXinv[0][1] = -aX[0][1] / determinantX; + m_aXinv[1][0] = -aX[1][0] / determinantX; + m_aXinv[1][1] = aX[0][0] / determinantX; + + double determinantY = aY[0][0] * aY[1][1] - aY[0][1] * aY[1][0]; + + if (determinantY == 0) { + error("Reco y matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!"); + return false; + } + + m_aYinv[0][0] = aY[1][1] / determinantY; + m_aYinv[0][1] = -aY[0][1] / determinantY; + m_aYinv[1][0] = -aY[1][0] / determinantY; + m_aYinv[1][1] = aY[0][0] / determinantY; + + return true; +} \ No newline at end of file diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.h b/src/algorithms/fardetectors/MatrixTransferStatic.h index b823e37e51..291ec1aa7f 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.h +++ b/src/algorithms/fardetectors/MatrixTransferStatic.h @@ -21,6 +21,7 @@ namespace eicrecon { using MatrixTransferStaticAlgorithm = algorithms::Algorithm< algorithms::Input< + edm4hep::MCParticleCollection, edm4hep::MCParticleCollection, edm4eic::TrackerHitCollection >, @@ -41,11 +42,24 @@ namespace eicrecon { "Apply matrix method reconstruction to hits."} {} void init() final; - void process(const Input&, const Output&) const final; + void process(const Input&, const Output&); + bool initalizeMatrix(const edm4hep::MCParticleCollection&); private: const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()}; const dd4hep::rec::CellIDPositionConverter* m_converter{algorithms::GeoSvc::instance().cellIDPositionConverter()}; + std::once_flag m_initBeamE; + bool m_initialized{false}; + double m_nomMomentum{0.0}; + double m_local_x_offset{0.0}; + double m_local_y_offset{0.0}; + double m_local_x_slope_offset{0.0}; + double m_local_y_slope_offset{0.0}; + + double m_aXinv[2][2] = {{0.0, 0.0}, + {0.0, 0.0}}; + double m_aYinv[2][2] = {{0.0, 0.0}, + {0.0, 0.0}}; }; } diff --git a/src/detectors/FOFFMTRK/FOFFMTRK.cc b/src/detectors/FOFFMTRK/FOFFMTRK.cc index 24a37ee54c..47c95a1697 100644 --- a/src/detectors/FOFFMTRK/FOFFMTRK.cc +++ b/src/detectors/FOFFMTRK/FOFFMTRK.cc @@ -67,7 +67,13 @@ void InitPlugin(JApplication *app) { recon_cfg.readout = "ForwardOffMTrackerRecHits"; - app->Add(new JOmniFactoryGeneratorT("ForwardOffMRecParticles",{"MCParticles","ForwardOffMTrackerRecHits"},{"ForwardOffMRecParticles"},recon_cfg,app)); + app->Add(new JOmniFactoryGeneratorT( + "ForwardOffMRecParticles", + {"MCBeamProtons","MCScatteredProtons","ForwardOffMTrackerRecHits"}, + {"ForwardOffMRecParticles"}, + recon_cfg, + app + )); } } diff --git a/src/detectors/RPOTS/RPOTS.cc b/src/detectors/RPOTS/RPOTS.cc index cf76ce7781..f3ab8fcb16 100644 --- a/src/detectors/RPOTS/RPOTS.cc +++ b/src/detectors/RPOTS/RPOTS.cc @@ -68,7 +68,13 @@ void InitPlugin(JApplication *app) { recon_cfg.readout = "ForwardRomanPotRecHits"; - app->Add(new JOmniFactoryGeneratorT("ForwardRomanPotRecParticles",{"MCParticles","ForwardRomanPotRecHits"},{"ForwardRomanPotRecParticles"},recon_cfg,app)); + app->Add(new JOmniFactoryGeneratorT( + "ForwardRomanPotRecParticles", + {"MCBeamProtons","MCScatteredProtons","ForwardRomanPotRecHits"}, + {"ForwardRomanPotRecParticles"}, + recon_cfg, + app + )); } } diff --git a/src/factories/fardetectors/MatrixTransferStatic_factory.h b/src/factories/fardetectors/MatrixTransferStatic_factory.h index f2bfd6943f..35506052df 100644 --- a/src/factories/fardetectors/MatrixTransferStatic_factory.h +++ b/src/factories/fardetectors/MatrixTransferStatic_factory.h @@ -25,7 +25,8 @@ class MatrixTransferStatic_factory : private: std::unique_ptr m_algo; - PodioInput m_mcparts_input {this}; + PodioInput m_mcbeamp_input {this}; + PodioInput m_mcscatp_input {this}; PodioInput m_hits_input {this}; PodioOutput m_tracks_output {this}; @@ -67,7 +68,7 @@ class MatrixTransferStatic_factory : } void Process(int64_t run_number, uint64_t event_number) { - m_algo->process({m_mcparts_input(), m_hits_input()}, {m_tracks_output().get()}); + m_algo->process({m_mcbeamp_input(), m_mcscatp_input(), m_hits_input()}, {m_tracks_output().get()}); } }; From cc9c8dcfd3af3fdb44b84907a98a7c147e2250f6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 13 Aug 2024 12:07:57 +0000 Subject: [PATCH 2/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/fardetectors/MatrixTransferStatic.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index 9319072e3a..092f9360be 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.cc +++ b/src/algorithms/fardetectors/MatrixTransferStatic.cc @@ -61,7 +61,7 @@ void eicrecon::MatrixTransferStatic::process( debug("No valid matrix can be identified for reconstruction"); return; } - + //---- begin Reconstruction code ---- edm4hep::Vector3f goodHit[2] = {{0.0,0.0,0.0},{0.0,0.0,0.0}}; @@ -231,7 +231,7 @@ bool eicrecon::MatrixTransferStatic::initalizeMatrix(const edm4hep::MCParticleCo m_local_y_offset = -0.00028343; m_local_x_slope_offset = -0.218525084; m_local_y_slope_offset = -0.00015321; - + m_nomMomentum = 100.0; } @@ -305,4 +305,4 @@ bool eicrecon::MatrixTransferStatic::initalizeMatrix(const edm4hep::MCParticleCo m_aYinv[1][1] = aY[0][0] / determinantY; return true; -} \ No newline at end of file +} From ff19849a03de378a73a60ba0be05c2b507bb10d3 Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 13 Aug 2024 13:13:35 +0100 Subject: [PATCH 3/6] Removed and moved variables --- .../fardetectors/MatrixTransferStatic.cc | 24 +++++-------------- 1 file changed, 6 insertions(+), 18 deletions(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index 092f9360be..0fad04581f 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.cc +++ b/src/algorithms/fardetectors/MatrixTransferStatic.cc @@ -22,6 +22,12 @@ void eicrecon::MatrixTransferStatic::init() { + m_nomMomentum = m_cfg.nomMomentum; //extract the nominal value first -- will be overwritten by MCParticle + m_local_x_offset = m_cfg.local_x_offset; + m_local_y_offset = m_cfg.local_y_offset; + m_local_x_slope_offset = m_cfg.local_x_slope_offset; + m_local_y_slope_offset = m_cfg.local_y_slope_offset; + } void eicrecon::MatrixTransferStatic::process( @@ -31,24 +37,6 @@ void eicrecon::MatrixTransferStatic::process( const auto [beamP,scatP,rechits] = input; auto [outputParticles] = output; - std::vector> aX(m_cfg.aX); - std::vector> aY(m_cfg.aY); - - //----- Define constants here ------ - double m_aXinv[2][2] = {{0.0, 0.0}, - {0.0, 0.0}}; - double m_aYinv[2][2] = {{0.0, 0.0}, - {0.0, 0.0}}; - - double nomMomentum = m_cfg.nomMomentum; //extract the nominal value first -- will be overwritten by MCParticle - double local_x_offset = m_cfg.local_x_offset; - double local_y_offset = m_cfg.local_y_offset; - double local_x_slope_offset = m_cfg.local_x_slope_offset; - double local_y_slope_offset = m_cfg.local_y_slope_offset; - - double numBeamProtons = 0; - double runningMomentum = 0.0; - //Set beam energy from first MCBeamElectron, using std::call_once std::call_once(m_initBeamE,[&](){ From 9d26e53b78675206c7fe49d268a83d68262d9be5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 13 Aug 2024 12:13:47 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/fardetectors/MatrixTransferStatic.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index 0fad04581f..7aa3fa99e0 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.cc +++ b/src/algorithms/fardetectors/MatrixTransferStatic.cc @@ -27,7 +27,7 @@ void eicrecon::MatrixTransferStatic::init() { m_local_y_offset = m_cfg.local_y_offset; m_local_x_slope_offset = m_cfg.local_x_slope_offset; m_local_y_slope_offset = m_cfg.local_y_slope_offset; - + } void eicrecon::MatrixTransferStatic::process( From ce9e4e446a33a406683fd436a111eb591c78ab3b Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 13 Aug 2024 13:28:50 +0100 Subject: [PATCH 5/6] Changed variables to have member tag --- src/algorithms/fardetectors/MatrixTransferStatic.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index 7aa3fa99e0..c4ab4f5e71 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.cc +++ b/src/algorithms/fardetectors/MatrixTransferStatic.cc @@ -106,8 +106,8 @@ void eicrecon::MatrixTransferStatic::process( // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit // trajectory - double XL[2] = {goodHit[0].x - local_x_offset, goodHit[1].x - local_x_offset}; - double YL[2] = {goodHit[0].y - local_y_offset, goodHit[1].y - local_y_offset}; + double XL[2] = {goodHit[0].x - m_local_x_offset, goodHit[1].x - m_local_x_offset}; + double YL[2] = {goodHit[0].y - m_local_y_offset, goodHit[1].y - m_local_y_offset}; double base = goodHit[1].z - goodHit[0].z; @@ -117,9 +117,9 @@ void eicrecon::MatrixTransferStatic::process( else{ double Xip[2] = {0.0, 0.0}; - double Xrp[2] = {XL[1], ((XL[1] - XL[0]) / (base))/dd4hep::mrad - local_x_slope_offset}; + double Xrp[2] = {XL[1], ((XL[1] - XL[0]) / (base))/dd4hep::mrad - m_local_x_slope_offset}; double Yip[2] = {0.0, 0.0}; - double Yrp[2] = {YL[1], ((YL[1] - YL[0]) / (base))/dd4hep::mrad - local_y_slope_offset}; + double Yrp[2] = {YL[1], ((YL[1] - YL[0]) / (base))/dd4hep::mrad - m_local_y_slope_offset}; // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the // Polar Angle and deltaP at the IP From feaa087ecee7582b84d5bf0b81f5f25f1358a3f7 Mon Sep 17 00:00:00 2001 From: simonge Date: Thu, 15 Aug 2024 16:24:21 +0100 Subject: [PATCH 6/6] Reverted where the proton filter is selected from --- .../fardetectors/MatrixTransferStatic.cc | 26 ++++++++++++++----- .../fardetectors/MatrixTransferStatic.h | 2 +- src/detectors/FOFFMTRK/FOFFMTRK.cc | 2 +- src/detectors/RPOTS/RPOTS.cc | 2 +- .../MatrixTransferStatic_factory.h | 5 ++-- 5 files changed, 25 insertions(+), 12 deletions(-) diff --git a/src/algorithms/fardetectors/MatrixTransferStatic.cc b/src/algorithms/fardetectors/MatrixTransferStatic.cc index c4ab4f5e71..2b6b2b9632 100644 --- 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 index 291ec1aa7f..58e8f4ed9c 100644 --- a/src/algorithms/fardetectors/MatrixTransferStatic.h +++ b/src/algorithms/fardetectors/MatrixTransferStatic.h @@ -13,6 +13,7 @@ #include #include #include +#include #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 index 47c95a1697..ab479c7ba7 100644 --- a/src/detectors/FOFFMTRK/FOFFMTRK.cc +++ b/src/detectors/FOFFMTRK/FOFFMTRK.cc @@ -69,7 +69,7 @@ void InitPlugin(JApplication *app) { app->Add(new JOmniFactoryGeneratorT( "ForwardOffMRecParticles", - {"MCBeamProtons","MCScatteredProtons","ForwardOffMTrackerRecHits"}, + {"MCParticles","ForwardOffMTrackerRecHits"}, {"ForwardOffMRecParticles"}, recon_cfg, app diff --git a/src/detectors/RPOTS/RPOTS.cc b/src/detectors/RPOTS/RPOTS.cc index f3ab8fcb16..4c643cb346 100644 --- a/src/detectors/RPOTS/RPOTS.cc +++ b/src/detectors/RPOTS/RPOTS.cc @@ -70,7 +70,7 @@ void InitPlugin(JApplication *app) { app->Add(new JOmniFactoryGeneratorT( "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 index 35506052df..f2bfd6943f 100644 --- 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 m_algo; - PodioInput m_mcbeamp_input {this}; - PodioInput m_mcscatp_input {this}; + PodioInput m_mcparts_input {this}; PodioInput m_hits_input {this}; PodioOutput 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()}); } };