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 <podio/ObjectID.h>
0013 #include <algorithm>
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 with reconstructed electron
0057   const auto ef_assoc = std::find_if(rcassoc->begin(), rcassoc->end(), [&ef_coll](const auto& a) {
0058     return a.getSim().getObjectID() == ef_coll[0].getObjectID();
0059   });
0060
0061   // Check to see if the associated reconstructed particle is available
0062   if (ef_assoc == rcassoc->end()) {
0063     trace("Truth scattered electron not in reconstructed particles");
0064     return;
0065   }
0066
0067   // Get the reconstructed electron object
0068   const auto ef_rc{(*ef_assoc).getRec()};
0069   const auto ef_rc_id{ef_rc.getObjectID()};
0070
0071   // Use these to compute the E-Pz
0072   // This is for development of the EMinusPz
0073   // algorithm
0074   PxPyPzMVector vScatteredElectron;
0075   PxPyPzMVector vHadronicFinalState;
0076   PxPyPzMVector vHadron;
0077
0078   // Loop over reconstructed particles to
0079   // get outgoing scattered electron.
0080   // Use the true scattered electron from the
0081   // MC information
0082   std::vector<PxPyPzEVector> electrons;
0083   for (const auto& p : *rcparts) {
0084     if (p.getObjectID() == ef_rc_id) {
0085
0086       output_electrons->push_back(p);
0087       static const auto m_electron = m_particleSvc.particle(11).mass;
0088       vScatteredElectron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0089                                         m_electron);
0090       electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z,
0091                              p.getEnergy());
0092       // break; NOTE: if we are not computing E-Pz
0093       // we can safely break here and save precious CPUs
0094     } else {
0095       // Compute the sum hadronic final state
0096       static const auto m_pion = m_particleSvc.particle(211).mass;
0097       vHadron.SetCoordinates(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_pion);
0098       vHadronicFinalState += vHadron;
0099     }
0100   }
0101
0102   // If no scattered electron was found, too bad
0103   if (electrons.empty()) {
0104     trace("No Truth scattered electron found");
0105     return;
0106   }
0107
0108   // Just for debug/development
0109   // report the computed E-Pz for the chosen electron
0110   double EPz = (vScatteredElectron + vHadronicFinalState).E() -
0111                (vScatteredElectron + vHadronicFinalState).Pz();
0112   trace("We found {} scattered electrons using Truth association", electrons.size());
0113   trace("TRUTH scattered electron has E-Pz = {}", EPz);
0114   trace("TRUTH scattered electron has Pxyz=( {}, {}, {} ) and E/p = {}", electrons[0].Px(),
0115         electrons[0].Py(), electrons[0].Pz(), (electrons[0].E() / electrons[0].P()));
0116 } // process
0117
0118 } // namespace eicrecon