algorithms/calorimetry/EnergyPositionClusterMerger.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten
0003
0004 #include "algorithms/calorimetry/EnergyPositionClusterMerger.h"
0005
0006 #include <edm4hep/Vector3f.h>
0007 #include <edm4hep/utils/vector_utils.h>
0008 #include <fmt/core.h>
0009 #include <podio/ObjectID.h>
0010 #include <cmath>
0011 #include <cstddef>
0012 #include <gsl/pointers>
0013 #include <limits>
0014 #include <vector>
0015
0016 #include "algorithms/calorimetry/EnergyPositionClusterMergerConfig.h"
0017
0018 namespace eicrecon {
0019
0020 void EnergyPositionClusterMerger::process(const Input& input, const Output& output) const {
0021
0022   const auto [energy_clus, energy_assoc, pos_clus, pos_assoc] = input;
0023   auto [merged_clus, merged_assoc]                            = output;
0024
0025   debug("Merging energy and position clusters for new event");
0026
0027   if (energy_clus->size() == 0 && pos_clus->size() == 0) {
0028     debug("Nothing to do for this event, returning...");
0029     return;
0030   }
0031
0032   std::vector<bool> consumed(energy_clus->size(), false);
0033
0034   // use position clusters as starting point
0035   for (const auto& pc : *pos_clus) {
0036
0037     trace(" --> Processing position cluster {}, energy: {}", pc.getObjectID().index,
0038           pc.getEnergy());
0039
0040     // check if we find a good match
0041     int best_match    = -1;
0042     double best_delta = std::numeric_limits<double>::max();
0043     for (std::size_t ie = 0; ie < energy_clus->size(); ++ie) {
0044       if (consumed[ie]) {
0045         continue;
0046       }
0047
0048       const auto& ec = (*energy_clus)[ie];
0049
0050       trace("  --> Evaluating energy cluster {}, energy: {}", ec.getObjectID().index,
0051             ec.getEnergy());
0052
0053       // 1. stop if not within tolerance
0054       //    (make sure to handle rollover of phi properly)
0055       const double de_rel = std::abs((pc.getEnergy() - ec.getEnergy()) / ec.getEnergy());
0056       const double deta =
0057           std::abs(edm4hep::utils::eta(pc.getPosition()) - edm4hep::utils::eta(ec.getPosition()));
0058       // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow
0059       // for phi rollovers
0060       const double dphi = edm4hep::utils::angleAzimuthal(pc.getPosition()) -
0061                           edm4hep::utils::angleAzimuthal(ec.getPosition());
0062       const double dsphi = std::abs(sin(0.5 * dphi));
0063       if ((m_cfg.energyRelTolerance > 0 && de_rel > m_cfg.energyRelTolerance) ||
0064           (m_cfg.etaTolerance > 0 && deta > m_cfg.etaTolerance) ||
0065           (m_cfg.phiTolerance > 0 && dsphi > sin(0.5 * m_cfg.phiTolerance))) {
0066         continue;
0067       }
0068       // --> if we get here, we have a match within tolerance. Now treat the case
0069       //     where we have multiple matches. In this case take the one with the closest
0070       //     energies.
0071       // 2. best match?
0072       const double delta = std::abs(pc.getEnergy() - ec.getEnergy());
0073       if (delta < best_delta) {
0074         best_delta = delta;
0075         best_match = ie;
0076       }
0077     }
0078
0079     // Create a merged cluster if we find a good match
0080     if (best_match >= 0) {
0081
0082       const auto& ec = (*energy_clus)[best_match];
0083
0084       auto new_clus = merged_clus->create();
0085       new_clus.setEnergy(ec.getEnergy());
0086       new_clus.setEnergyError(ec.getEnergyError());
0087       new_clus.setTime(pc.getTime());
0088       new_clus.setNhits(pc.getNhits() + ec.getNhits());
0089       new_clus.setPosition(pc.getPosition());
0090       new_clus.setPositionError(pc.getPositionError());
0091       new_clus.addToClusters(pc);
0092       new_clus.addToClusters(ec);
0093
0094       trace("   --> Found matching energy cluster {}, energy: {}", ec.getObjectID().index,
0095             ec.getEnergy());
0096       trace("   --> Created a new combined cluster {}, energy: {}", new_clus.getObjectID().index,
0097             new_clus.getEnergy());
0098
0099       // find association from energy cluster
0100       auto ea = energy_assoc->begin();
0101       for (; ea != energy_assoc->end(); ++ea) {
0102         if (ea->getRec() == ec) {
0103           break;
0104         }
0105       }
0106       // find association from position cluster if different
0107       auto pa = pos_assoc->begin();
0108       for (; pa != pos_assoc->end(); ++pa) {
0109         if (pa->getRec() == pc) {
0110           break;
0111         }
0112       }
0113       if (ea != energy_assoc->end() || pa != pos_assoc->end()) {
0114         // we must write an association
0115         if (ea != energy_assoc->end() && pa != pos_assoc->end()) {
0116           // we have two associations
0117           if (pa->getSimID() == ea->getSimID()) {
0118             // both associations agree on the MCParticles entry
0119             auto clusterassoc = merged_assoc->create();
0120             clusterassoc.setRecID(new_clus.getObjectID().index);
0121             clusterassoc.setSimID(ea->getSimID());
0122             clusterassoc.setWeight(1.0);
0123             clusterassoc.setRec(new_clus);
0124             clusterassoc.setSim(ea->getSim());
0125           } else {
0126             // both associations disagree on the MCParticles entry
0127             debug("   --> Two associations added to {} and {}", ea->getSimID(), pa->getSimID());
0128             auto clusterassoc1 = merged_assoc->create();
0129             clusterassoc1.setRecID(new_clus.getObjectID().index);
0130             clusterassoc1.setSimID(ea->getSimID());
0131             clusterassoc1.setWeight(0.5);
0132             clusterassoc1.setRec(new_clus);
0133             clusterassoc1.setSim(ea->getSim());
0134             auto clusterassoc2 = merged_assoc->create();
0135             clusterassoc2.setRecID(new_clus.getObjectID().index);
0136             clusterassoc2.setSimID(pa->getSimID());
0137             clusterassoc2.setWeight(0.5);
0138             clusterassoc2.setRec(new_clus);
0139             clusterassoc2.setSim(pa->getSim());
0140           }
0141         } else if (ea != energy_assoc->end()) {
0142           // no position association
0143           debug("   --> Only added energy cluster association to {}", ea->getSimID());
0144           auto clusterassoc = merged_assoc->create();
0145           clusterassoc.setRecID(new_clus.getObjectID().index);
0146           clusterassoc.setSimID(ea->getSimID());
0147           clusterassoc.setWeight(1.0);
0148           clusterassoc.setRec(new_clus);
0149           clusterassoc.setSim(ea->getSim());
0150         } else if (pa != pos_assoc->end()) {
0151           // no energy association
0152           debug("   --> Only added position cluster association to {}", pa->getSimID());
0153           auto clusterassoc = merged_assoc->create();
0154           clusterassoc.setRecID(new_clus.getObjectID().index);
0155           clusterassoc.setSimID(pa->getSimID());
0156           clusterassoc.setWeight(1.0);
0157           clusterassoc.setRec(new_clus);
0158           clusterassoc.setSim(pa->getSim());
0159         }
0160       }
0161
0162       // label our energy cluster as consumed
0163       consumed[best_match] = true;
0164
0165       debug("  Matched position cluster {} with energy cluster {}", pc.getObjectID().index,
0166             ec.getObjectID().index);
0167       debug("  - Position cluster: (E: {}, phi: {}, z: {})", pc.getEnergy(),
0168             edm4hep::utils::angleAzimuthal(pc.getPosition()), pc.getPosition().z);
0169       debug("  - Energy cluster: (E: {}, phi: {}, z: {})", ec.getEnergy(),
0170             edm4hep::utils::angleAzimuthal(ec.getPosition()), ec.getPosition().z);
0171       debug("  ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.getEnergy(),
0172             edm4hep::utils::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z);
0173
0174     } else {
0175
0176       debug("  Unmatched position cluster {}", pc.getObjectID().index);
0177     }
0178   }
0179 }
0180
0181 } // namespace eicrecon