algorithms/reco/ElectronReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Daniel Brandenburg
0003 #include "ElectronReconstruction.h"
0004
0005 #include <edm4eic/ClusterCollection.h>
0006 #include <edm4eic/ReconstructedParticleCollection.h>
0007 #include <edm4hep/utils/vector_utils.h>
0008 #include <fmt/core.h>
0009 #include <podio/RelationRange.h>
0010
0011 #include "algorithms/reco/ElectronReconstructionConfig.h"
0012
0013 namespace eicrecon {
0014
0015     void ElectronReconstruction::init(std::shared_ptr<spdlog::logger> logger) {
0016         m_log = logger;
0017     }
0018
0019     std::unique_ptr<edm4eic::ReconstructedParticleCollection> ElectronReconstruction::execute(
0020             const edm4eic::ReconstructedParticleCollection *rcparts
0021             ) {
0022
0023         // Some obvious improvements:
0024         // - E/p cut from real study optimized for electron finding and hadron rejection
0025         // - use of any HCAL info?
0026         // - check for duplicates?
0027
0028         // output container
0029         auto out_electrons = std::make_unique<edm4eic::ReconstructedParticleCollection>();
0030         out_electrons->setSubsetCollection(); // out_electrons is a subset of the ReconstructedParticles collection
0031
0032         for (const auto particle : *rcparts) {
0033             // if we found a reco particle then test for electron compatibility
0034             if (particle.getClusters().size() == 0) {
0035                 continue;
0036             }
0037             if (particle.getCharge() == 0) { // Skip over photons/other particles without a track
0038                 continue;
0039             }
0040             double E = particle.getClusters()[0].getEnergy();
0041             double p = edm4hep::utils::magnitude(particle.getMomentum());
0042             double EOverP = E / p;
0043
0044             m_log->trace("ReconstructedElectron: Energy={} GeV, p={} GeV, E/p = {} for PDG (from truth): {}", E, p, EOverP, particle.getPDG());
0045             // Apply the E/p cut here to select electons
0046             if (EOverP >= m_cfg.min_energy_over_momentum && EOverP <= m_cfg.max_energy_over_momentum) {
0047                 out_electrons->push_back(particle);
0048             }
0049
0050         }
0051     m_log->debug("Found {} electron candidates", out_electrons->size());
0052     return out_electrons;
0053     }
0054 }