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(
0037       const ScatteredElectronsEMinusPz::Input& input,
0038       const ScatteredElectronsEMinusPz::Output& output
0039   ) const {
0040     auto [rcparts, rcele] = input;
0041     auto [out_electrons] = output;
0042
0043     static const auto m_electron = m_particleSvc.particle(11).mass;
0044     static const auto m_pion = m_particleSvc.particle(211).mass;
0045
0046     // this map will store intermediate results
0047     // so that we can sort them before filling output
0048     // collection
0049     std::map<double, edm4eic::ReconstructedParticle, std::greater<double>> scatteredElectronsMap;
0050
0051     out_electrons->setSubsetCollection();
0052
0053     trace("We have {} candidate electrons", rcele->size());
0054
0055     // Lorentz Vector for the scattered electron,
0056     // hadronic final state, and individual hadron
0057     // We do it here to avoid creating objects inside the loops
0058     PxPyPzMVector  vScatteredElectron;
0059     PxPyPzMVector  vHadronicFinalState;
0060     PxPyPzMVector  vHadron;
0061
0062     for ( const auto& e: *rcele ) {
0063       // Do not cut on charge to account for charge-symmetric background
0064
0065       // reset the HadronicFinalState
0066       vHadronicFinalState.SetCoordinates(0, 0, 0, 0);
0067
0068       // Set a vector for the electron we are considering now
0069       vScatteredElectron.SetCoordinates(
0070           e.getMomentum().x,
0071           e.getMomentum().y,
0072           e.getMomentum().z,
0073           m_electron
0074         );
0075
0076       // Loop over reconstructed particles to
0077       // sum hadronic final state
0078       for (const auto& p: *rcparts) {
0079         // What we want is to add all reconstructed particles
0080         // except the one we are currently considering as the
0081         // (scattered) electron candidate.
0082         if (p.getObjectID() != e.getObjectID()) {
0083           vHadron.SetCoordinates(
0084               p.getMomentum().x,
0085               p.getMomentum().y,
0086               p.getMomentum().z,
0087               m_pion // Assume pion for hadronic state
0088             );
0089
0090           // Sum hadronic final state
0091           vHadronicFinalState += vHadron;
0092         } else {
0093           trace("Skipping electron in hadronic final state");
0094         }
0095       } // hadron loop (reconstructed particles)
0096
0097       // Calculate the E-Pz for this electron
0098       // + hadron final state combination
0099       // For now we keep all electron
0100       // candidates but we will rank them by their E-Pz
0101       double EPz=(vScatteredElectron+vHadronicFinalState).E()
0102               - (vScatteredElectron+vHadronicFinalState).Pz();
0103       trace("\tE-Pz={}", EPz);
0104       trace("\tScatteredElectron has Pxyz=( {}, {}, {} )", e.getMomentum().x, e.getMomentum().y, e.getMomentum().z);
0105
0106       // Store the result of this calculation
0107       scatteredElectronsMap[ EPz ] = e;
0108     } // electron loop
0109
0110     trace("Selecting candidates with {} < E-Pz < {}", m_cfg.minEMinusPz, m_cfg.maxEMinusPz);
0111
0112     bool first = true;
0113     // map defined with std::greater<> will be iterated in descending order by the key
0114     for (auto kv : scatteredElectronsMap) {
0115
0116       double EMinusPz = kv.first;
0117       // Do not save electron candidates that
0118       // are not within range
0119       if ( EMinusPz > m_cfg.maxEMinusPz
0120         || EMinusPz < m_cfg.minEMinusPz ){
0121         continue;
0122       }
0123
0124       // For logging and development
0125       // report the highest E-Pz candidate chosen
0126       if ( first ){
0127         trace("Max E-Pz Candidate:");
0128         trace("\tE-Pz={}", EMinusPz);
0129         trace("\tScatteredElectron has Pxyz=( {}, {}, {} )", kv.second.getMomentum().x, kv.second.getMomentum().y, kv.second.getMomentum().z);
0130         first = false;
0131       }
0132       out_electrons->push_back( kv.second );
0133     } // reverse loop on scatteredElectronsMap
0134   }
0135
0136 } // namespace eicrecon