calorimetry/todo/EnergyPositionClusterMerger.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten
0003
0004 #include <limits>
0005 #include <numbers>
0006
0007 #include <fmt/format.h>
0008 // Gaudi
0009 #include "Gaudi/Property.h"
0010 #include "GaudiAlg/GaudiAlgorithm.h"
0011 #include "GaudiAlg/GaudiTool.h"
0012 #include "GaudiAlg/Transformer.h"
0013 #include "GaudiKernel/PhysicalConstants.h"
0014
0015 #include "JugBase/DataHandle.h"
0016
0017 // Event Model related classes
0018 #include "edm4eic/ClusterCollection.h"
0019 #include "edm4hep/utils/vector_utils.h"
0020
0021 using namespace Gaudi::Units;
0022
0023 namespace Jug::Reco {
0024
0025 /** Simple algorithm to merge the energy measurement from cluster1 with the position
0026  * measurement of cluster2 (in case matching clusters are found). If not, it will
0027  * propagate the raw cluster from cluster1 or cluster2
0028  *
0029  * Matching occurs based on the cluster phi, z and E variables, with tolerances
0030  * defined in the options file. A negative tolerance effectively disables
0031  * a check. The energy tolerance is defined as a relative number (e.g. .1)
0032  *
0033  * In case of ambiguity the closest cluster is merged.
0034  *
0035  * \ingroup reco
0036  */
0037 class EnergyPositionClusterMerger : public GaudiAlgorithm {
0038 private:
0039   // Input
0040   DataHandle<edm4eic::ClusterCollection> m_energyClusters{"energyClusters", Gaudi::DataHandle::Reader, this};
0041   DataHandle<edm4eic::ClusterCollection> m_positionClusters{"positionClusters", Gaudi::DataHandle::Reader, this};
0042   // Output
0043   DataHandle<edm4eic::ClusterCollection> m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this};
0044   // Negative values mean the tolerance check is disabled
0045   Gaudi::Property<double> m_zToleranceUnits{this, "zTolerance", -1 * cm};
0046   Gaudi::Property<double> m_phiToleranceUnits{this, "phiTolerance", 20 * degree};
0047   Gaudi::Property<double> m_energyRelTolerance{this, "energyRelTolerance", 0.3};
0048   // Unitless (GeV/mm/ns/rad) versions of these tolerances
0049   double m_zTolerance{0};
0050   double m_phiTolerance{0};
0051
0052 public:
0053   EnergyPositionClusterMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0054     declareProperty("energyClusters", m_energyClusters, "Cluster collection with good energy precision");
0055     declareProperty("positionClusters", m_positionClusters, "Cluster collection with good position precision");
0056     declareProperty("outputClusters", m_outputClusters, "");
0057   }
0058
0059   StatusCode initialize() override {
0060     m_zTolerance   = m_zToleranceUnits / mm;
0061     m_phiTolerance = m_phiTolerance / rad;
0062     return StatusCode::SUCCESS;
0063   }
0064
0065   StatusCode execute() override {
0066     // input
0067     const auto& e_clus   = *(m_energyClusters.get());
0068     const auto& pos_clus = *(m_positionClusters.get());
0069     // output
0070     auto& merged = *(m_outputClusters.createAndPut());
0071
0072     std::vector<bool> consumed(e_clus.size(), false);
0073
0074     // use position clusters as starting point
0075     for (const auto& pc : pos_clus) {
0076       // check if we find a good match
0077       int best_match    = -1;
0078       double best_delta = std::numeric_limits<double>::max();
0079       for (size_t ie = 0; ie < e_clus.size(); ++ie) {
0080         if (consumed[ie]) {
0081           continue;
0082         }
0083         const auto& ec = e_clus[ie];
0084         // 1. stop if not within tolerance
0085         //    (make sure to handle rollover of phi properly)
0086         double dphi = edm4hep::utils::angleAzimuthal(pc.getPosition()) - edm4hep::utils::angleAzimuthal(ec.getPosition());
0087         if (std::abs(dphi) > M_PI) {
0088           dphi = std::abs(dphi) - M_PI;
0089         }
0090         if ((m_energyRelTolerance > 0 &&
0091              (ec.getEnergy() <= 0 || fabs((pc.getEnergy() - ec.getEnergy()) / ec.getEnergy()) > m_energyRelTolerance)) ||
0092             (m_zTolerance > 0 && fabs(pc.getPosition().z - ec.getPosition().z) > m_zTolerance) ||
0093             (m_phiTolerance > 0 && dphi > m_phiTolerance)) {
0094           continue;
0095         }
0096         // --> if we get here, we have a match within tolerance. Now treat the case
0097         //     where we have multiple matches. In this case take the one with the closest
0098         //     energies.
0099         // 2. best match?
0100         const double delta = fabs(pc.getEnergy() - ec.getEnergy());
0101         if (delta < best_delta) {
0102           best_delta = delta;
0103           best_match = ie;
0104         }
0105       }
0106       // Create a merged cluster if we find a good match
0107       if (best_match >= 0) {
0108         const auto& ec = e_clus[best_match];
0109         auto new_clus  = merged.create();
0110         new_clus.setEnergy(ec.getEnergy());
0111         new_clus.setEnergyError(ec.getEnergyError());
0112         new_clus.setTime(pc.getTime());
0113         new_clus.setNhits(pc.getNhits() + ec.getNhits());
0114         new_clus.setPosition(pc.getPosition());
0115         new_clus.setPositionError(pc.getPositionError());
0116         new_clus.addToClusters(pc);
0117         new_clus.addToClusters(ec);
0118         // label our energy cluster as consumed
0119         consumed[best_match] = true;
0120         if (msgLevel(MSG::DEBUG)) {
0121           debug() << fmt::format("Matched position cluster {} with energy cluster {}\n", pc.id(), ec.id()) << endmsg;
0122           debug() << fmt::format("  - Position cluster: (E: {}, phi: {}, z: {})", pc.getEnergy(),
0123                                  edm4hep::utils::angleAzimuthal(pc.getPosition()), pc.getPosition().z)
0124                   << endmsg;
0125           debug() << fmt::format("  - Energy cluster: (E: {}, phi: {}, z: {})", ec.getEnergy(),
0126                                  edm4hep::utils::angleAzimuthal(ec.getPosition()), ec.getPosition().z)
0127                   << endmsg;
0128           debug() << fmt::format("  ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.getEnergy(),
0129                                  edm4hep::utils::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z)
0130                   << endmsg;
0131         }
0132       }
0133     }
0134     // That's all!
0135
0136     return StatusCode::SUCCESS;
0137   }
0138 };
0139 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0140 DECLARE_COMPONENT(EnergyPositionClusterMerger)
0141
0142 } // namespace Jug::Reco