algorithms/reco/ScatteredElectronsTruth.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 Daniel Brandenburg
0003
0004 #include <Math/GenVector/LorentzVector.h>
0005 #include <Math/GenVector/PxPyPzE4D.h>
0006 #include <Math/GenVector/PxPyPzM4D.h>
0007 #include <Math/Vector4Dfwd.h>
0008 #include <edm4eic/MCRecoParticleAssociationCollection.h>
0009 #include <edm4eic/ReconstructedParticleCollection.h>
0010 #include <edm4hep/MCParticleCollection.h>
0011 #include <edm4hep/Vector3f.h>
0012 #include <fmt/core.h>
0013 #include <podio/ObjectID.h>
0014 #include <gsl/pointers>
0015 #include <vector>
0016
0017 #include "algorithms/reco/Beam.h"
0018 #include "algorithms/reco/ScatteredElectronsTruth.h"
0019
0020 using ROOT::Math::PxPyPzEVector;
0021 using ROOT::Math::PxPyPzMVector;
0022
0023 namespace eicrecon {
0024
0025 /**
0026    * @brief Initialize ScatteredElectronsTruth algorithm
0027    */
0028 void ScatteredElectronsTruth::init() {}
0029
0030 /**
0031    * @brief Selects the scattered electron based
0032    * on TRUTH information (Mc particle info).
0033    * Returns a collection of reconstructed particles
0034    * corresponding to the chosen Mc electrons
0035    *
0036    * @param input - McParticleCollection,
0037    *                ReconstructedParticleCollection,
0038    *                MCRecoParticleAssociationCollection
0039    * @param output - ReconstructedParticleCollection
0040    */
0041 void ScatteredElectronsTruth::process(const ScatteredElectronsTruth::Input& input,
0042                                       const ScatteredElectronsTruth::Output& output) const {
0043
0044   // get our input and outputs
0045   const auto [mcparts, rcparts, rcassoc] = input;
0046   auto [output_electrons]                = output;
0047   output_electrons->setSubsetCollection();
0048
0049   // Get first scattered electron
0050   const auto ef_coll = find_first_scattered_electron(mcparts);
0051   if (ef_coll.empty()) {
0052     trace("No truth scattered electron found");
0053     return;
0054   }
0055
0056   // Associate first scattered electron
0057   // with reconstructed electron
0058   auto ef_assoc = rcassoc->begin();
0059   for (; ef_assoc != rcassoc->end(); ++ef_assoc) {
0060     if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) {
0061       break;
0062     }
0063   }
0064
0065   // Check to see if the associated reconstructed
0066   // particle is available
0067   if (!(ef_assoc != rcassoc->end())) {
0068     trace("Truth scattered electron not in reconstructed particles");
0069     return;
0070   }
0071
0072   // Get the reconstructed electron object
0073   const auto ef_rc{ef_assoc->getRec()};
0074   const auto ef_rc_id{ef_rc.getObjectID()};
0075
0076   // Use these to compute the E-Pz
0077   // This is for development of the EMinusPz
0078   // algorithm
0079   PxPyPzMVector vScatteredElectron;
0080   PxPyPzMVector vHadronicFinalState;
0081   PxPyPzMVector vHadron;
0082
0083   // Loop over reconstructed particles to
0084   // get outgoing scattered electron.
0085   // Use the true scattered electron from the
0086   // MC information
0087   std::vector<PxPyPzEVector> electrons;
0088   for (const auto& p : *rcparts) {
0089     if (p.getObjectID() == ef_rc_id) {
0090
0091       output_electrons->push_back(p);
0092       static const auto m_electron = m_particleSvc.particle(11).mass;
0093       vScatteredElectron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0094                                         m_electron);
0095       electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0096                              p.getEnergy());
0097       // break; NOTE: if we are not computing E-Pz
0098       // we can safely break here and save precious CPUs
0099     } else {
0100       // Compute the sum hadronic final state
0101       static const auto m_pion = m_particleSvc.particle(211).mass;
0102       vHadron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_pion);
0103       vHadronicFinalState += vHadron;
0104     }
0105   }
0106
0107   // If no scattered electron was found, too bad
0108   if (electrons.empty()) {
0109     trace("No Truth scattered electron found");
0110     return;
0111   }
0112
0113   // Just for debug/development
0114   // report the computed E-Pz for the chosen electron
0115   double EPz = (vScatteredElectron + vHadronicFinalState).E() -
0116                (vScatteredElectron + vHadronicFinalState).Pz();
0117   trace("We found {} scattered electrons using Truth association", electrons.size());
0118   trace("TRUTH scattered electron has E-Pz = {}", EPz);
0119   trace("TRUTH scattered electron has Pxyz=( {}, {}, {} ) and E/p = {}", electrons[0].Px(),
0120         electrons[0].Py(), electrons[0].Pz(), (electrons[0].E() / electrons[0].P()));
0121 } // process
0122
0123 } // namespace eicrecon