File indexing completed on 2025-07-11 07:53:36
0001
0002
0003
0004 #include <Math/GenVector/LorentzVector.h>
0005 #include <Math/GenVector/PxPyPzM4D.h>
0006 #include <Math/Vector4Dfwd.h>
0007 #include <edm4eic/ReconstructedParticleCollection.h>
0008 #include <edm4hep/Vector3f.h>
0009 #include <fmt/core.h>
0010 #include <podio/ObjectID.h>
0011 #include <functional>
0012 #include <gsl/pointers>
0013 #include <map>
0014 #include <utility>
0015
0016 #include "algorithms/reco/ScatteredElectronsEMinusPz.h"
0017 #include "algorithms/reco/ScatteredElectronsEMinusPzConfig.h"
0018
0019 using ROOT::Math::PxPyPzEVector;
0020 using ROOT::Math::PxPyPzMVector;
0021
0022 namespace eicrecon {
0023
0024
0025
0026
0027 void ScatteredElectronsEMinusPz::init() {}
0028
0029
0030
0031
0032
0033
0034
0035
0036 void ScatteredElectronsEMinusPz::process(const ScatteredElectronsEMinusPz::Input& input,
0037 const ScatteredElectronsEMinusPz::Output& output) const {
0038 auto [rcparts, rcele] = input;
0039 auto [out_electrons] = output;
0040
0041 static const auto m_electron = m_particleSvc.particle(11).mass;
0042 static const auto m_pion = m_particleSvc.particle(211).mass;
0043
0044
0045
0046
0047 std::map<double, edm4eic::ReconstructedParticle, std::greater<>> scatteredElectronsMap;
0048
0049 out_electrons->setSubsetCollection();
0050
0051 trace("We have {} candidate electrons", rcele->size());
0052
0053
0054
0055
0056 PxPyPzMVector vScatteredElectron;
0057 PxPyPzMVector vHadronicFinalState;
0058 PxPyPzMVector vHadron;
0059
0060 for (const auto& e : *rcele) {
0061
0062
0063
0064 vHadronicFinalState.SetCoordinates(0, 0, 0, 0);
0065
0066
0067 vScatteredElectron.SetCoordinates(e.getMomentum().x, e.getMomentum().y, e.getMomentum().z,
0068 m_electron);
0069
0070
0071
0072 for (const auto& p : *rcparts) {
0073
0074
0075
0076 if (p.getObjectID() != e.getObjectID()) {
0077 vHadron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0078 m_pion
0079 );
0080
0081
0082 vHadronicFinalState += vHadron;
0083 } else {
0084 trace("Skipping electron in hadronic final state");
0085 }
0086 }
0087
0088
0089
0090
0091
0092 double EPz = (vScatteredElectron + vHadronicFinalState).E() -
0093 (vScatteredElectron + vHadronicFinalState).Pz();
0094 trace("\tE-Pz={}", EPz);
0095 trace("\tScatteredElectron has Pxyz=( {}, {}, {} )", e.getMomentum().x, e.getMomentum().y,
0096 e.getMomentum().z);
0097
0098
0099 scatteredElectronsMap[EPz] = e;
0100 }
0101
0102 trace("Selecting candidates with {} < E-Pz < {}", m_cfg.minEMinusPz, m_cfg.maxEMinusPz);
0103
0104 bool first = true;
0105
0106 for (auto kv : scatteredElectronsMap) {
0107
0108 double EMinusPz = kv.first;
0109
0110
0111 if (EMinusPz > m_cfg.maxEMinusPz || EMinusPz < m_cfg.minEMinusPz) {
0112 continue;
0113 }
0114
0115
0116
0117 if (first) {
0118 trace("Max E-Pz Candidate:");
0119 trace("\tE-Pz={}", EMinusPz);
0120 trace("\tScatteredElectron has Pxyz=( {}, {}, {} )", kv.second.getMomentum().x,
0121 kv.second.getMomentum().y, kv.second.getMomentum().z);
0122 first = false;
0123 }
0124 out_electrons->push_back(kv.second);
0125 }
0126 }
0127
0128 }