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(
0042       const ScatteredElectronsTruth::Input& input,
0043       const ScatteredElectronsTruth::Output& output) const {
0044
0045     // get our input and outputs
0046     const auto [mcparts, rcparts, rcassoc] = input;
0047     auto [output_electrons] = output;
0048     output_electrons->setSubsetCollection();
0049
0050
0051     // Get first scattered electron
0052     const auto ef_coll = find_first_scattered_electron(mcparts);
0053     if (ef_coll.size() == 0) {
0054       trace("No truth scattered electron found");
0055       return;
0056     }
0057
0058     // Associate first scattered electron
0059     // with reconstructed electron
0060     auto ef_assoc = rcassoc->begin();
0061     for (; ef_assoc != rcassoc->end(); ++ef_assoc) {
0062       if (ef_assoc->getSim().getObjectID() ==  ef_coll[0].getObjectID()) {
0063         break;
0064       }
0065     }
0066
0067     // Check to see if the associated reconstructed
0068     // particle is available
0069     if (!(ef_assoc != rcassoc->end())) {
0070       trace("Truth scattered electron not in reconstructed particles");
0071       return;
0072     }
0073
0074     // Get the reconstructed electron object
0075     const auto ef_rc{ef_assoc->getRec()};
0076     const auto ef_rc_id{ef_rc.getObjectID()};
0077
0078     // Use these to compute the E-Pz
0079     // This is for development of the EMinusPz
0080     // algorithm
0081     PxPyPzMVector  vScatteredElectron;
0082     PxPyPzMVector  vHadronicFinalState;
0083     PxPyPzMVector  vHadron;
0084
0085     // Loop over reconstructed particles to
0086     // get outgoing scattered electron.
0087     // Use the true scattered electron from the
0088     // MC information
0089     std::vector<PxPyPzEVector> electrons;
0090     for (const auto& p: *rcparts) {
0091       if (p.getObjectID() == ef_rc_id) {
0092
0093         output_electrons->push_back( p );
0094         static const auto m_electron = m_particleSvc.particle(11).mass;
0095         vScatteredElectron.SetCoordinates( p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, m_electron );
0096         electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, 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.size() == 0) {
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(
0120         "TRUTH scattered electron has Pxyz=( {}, {}, {} ) and E/p = {}",
0121         electrons[0].Px(),
0122         electrons[0].Py(),
0123         electrons[0].Pz(),
0124         (electrons[0].E()/electrons[0].P())
0125       );
0126   } // process
0127
0128 } // namespace eicrecon