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