algorithms/reco/ScatteredElectronsEMinusPz.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 Daniel Brandenburg, Dmitry Kalinkin
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    * @brief Initialize the ScatteredElectronsEMinusPz Algorithm
0026    */
0027 void ScatteredElectronsEMinusPz::init() {}
0028
0029 /**
0030    * @brief Produce a list of scattered electron candidates
0031    *
0032    * @param rcparts - input collection of all reconstructed particles
0033    * @param rcele  - input collection of all electron candidates
0034    * @return std::unique_ptr<edm4eic::ReconstructedParticleCollection>
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   // this map will store intermediate results
0045   // so that we can sort them before filling output
0046   // collection
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   // Lorentz Vector for the scattered electron,
0054   // hadronic final state, and individual hadron
0055   // We do it here to avoid creating objects inside the loops
0056   PxPyPzMVector vScatteredElectron;
0057   PxPyPzMVector vHadronicFinalState;
0058   PxPyPzMVector vHadron;
0059
0060   for (const auto& e : *rcele) {
0061     // Do not cut on charge to account for charge-symmetric background
0062
0063     // reset the HadronicFinalState
0064     vHadronicFinalState.SetCoordinates(0, 0, 0, 0);
0065
0066     // Set a vector for the electron we are considering now
0067     vScatteredElectron.SetCoordinates(e.getMomentum().x, e.getMomentum().y, e.getMomentum().z,
0068                                       m_electron);
0069
0070     // Loop over reconstructed particles to
0071     // sum hadronic final state
0072     for (const auto& p : *rcparts) {
0073       // What we want is to add all reconstructed particles
0074       // except the one we are currently considering as the
0075       // (scattered) electron candidate.
0076       if (p.getObjectID() != e.getObjectID()) {
0077         vHadron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0078                                m_pion // Assume pion for hadronic state
0079         );
0080
0081         // Sum hadronic final state
0082         vHadronicFinalState += vHadron;
0083       } else {
0084         trace("Skipping electron in hadronic final state");
0085       }
0086     } // hadron loop (reconstructed particles)
0087
0088     // Calculate the E-Pz for this electron
0089     // + hadron final state combination
0090     // For now we keep all electron
0091     // candidates but we will rank them by their E-Pz
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     // Store the result of this calculation
0099     scatteredElectronsMap[EPz] = e;
0100   } // electron loop
0101
0102   trace("Selecting candidates with {} < E-Pz < {}", m_cfg.minEMinusPz, m_cfg.maxEMinusPz);
0103
0104   bool first = true;
0105   // map defined with std::greater<> will be iterated in descending order by the key
0106   for (auto kv : scatteredElectronsMap) {
0107
0108     double EMinusPz = kv.first;
0109     // Do not save electron candidates that
0110     // are not within range
0111     if (EMinusPz > m_cfg.maxEMinusPz || EMinusPz < m_cfg.minEMinusPz) {
0112       continue;
0113     }
0114
0115     // For logging and development
0116     // report the highest E-Pz candidate chosen
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   } // reverse loop on scatteredElectronsMap
0126 }
0127
0128 } // namespace eicrecon