algorithms/calorimetry/EnergyPositionClusterMerger.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten
0003
0004 #pragma once
0005
0006 #include <limits>
0007
0008 #include <algorithms/algorithm.h>
0009 #include <fmt/format.h>
0010 #include <spdlog/spdlog.h>
0011 #include <edm4eic/ClusterCollection.h>
0012 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0013 #include <edm4hep/utils/vector_utils.h>
0014
0015 #include "algorithms/interfaces/WithPodConfig.h"
0016 #include "EnergyPositionClusterMergerConfig.h"
0017
0018 namespace eicrecon {
0019
0020   using EnergyPositionClusterMergerAlgorithm = algorithms::Algorithm<
0021     algorithms::Input<
0022       edm4eic::ClusterCollection,
0023       edm4eic::MCRecoClusterParticleAssociationCollection,
0024       edm4eic::ClusterCollection,
0025       edm4eic::MCRecoClusterParticleAssociationCollection
0026     >,
0027     algorithms::Output<
0028       edm4eic::ClusterCollection,
0029       edm4eic::MCRecoClusterParticleAssociationCollection
0030     >
0031   >;
0032
0033   /** Simple algorithm to merge the energy measurement from cluster1 with the position
0034   * measurement of cluster2 (in case matching clusters are found). If not, it will
0035   * propagate the raw cluster from cluster1 or cluster2
0036   *
0037   * Matching occurs based on the cluster phi, eta and E variables, with tolerances
0038   * defined in the options file. A negative tolerance effectively disables
0039   * a check. The energy tolerance is defined as a relative number (e.g. 0.1)
0040   *
0041   * In case of ambiguity the closest cluster is merged.
0042   *
0043   * \ingroup reco
0044   */
0045   class EnergyPositionClusterMerger
0046       : public EnergyPositionClusterMergerAlgorithm,
0047         public WithPodConfig<EnergyPositionClusterMergerConfig> {
0048
0049   public:
0050     EnergyPositionClusterMerger(std::string_view name)
0051       : EnergyPositionClusterMergerAlgorithm{name,
0052                             {"energyClusterCollection", "energyClusterAssociations",
0053                              "positionClusterCollection", "positionClusterAssociations"},
0054                             {"outputClusterCollection", "outputClusterAssociations"},
0055                             "Merge energy and position clusters if matching."} {}
0056
0057   public:
0058
0059     void init() { }
0060
0061     void process(const Input& input, const Output& output) const final {
0062
0063         const auto [energy_clus, energy_assoc, pos_clus, pos_assoc] = input;
0064         auto [merged_clus, merged_assoc] = output;
0065
0066         debug( "Merging energy and position clusters for new event" );
0067
0068         if (energy_clus->size() == 0 && pos_clus->size() == 0) {
0069             debug( "Nothing to do for this event, returning..." );
0070             return;
0071         }
0072
0073         std::vector<bool> consumed(energy_clus->size(), false);
0074
0075         // use position clusters as starting point
0076         for (const auto& pc : *pos_clus) {
0077
0078             trace(" --> Processing position cluster {}, energy: {}", pc.getObjectID().index, pc.getEnergy());
0079
0080             // check if we find a good match
0081             int best_match    = -1;
0082             double best_delta = std::numeric_limits<double>::max();
0083             for (size_t ie = 0; ie < energy_clus->size(); ++ie) {
0084                 if (consumed[ie]) {
0085                     continue;
0086                 }
0087
0088                 const auto& ec = (*energy_clus)[ie];
0089
0090                 trace("  --> Evaluating energy cluster {}, energy: {}", ec.getObjectID().index, ec.getEnergy());
0091
0092                 // 1. stop if not within tolerance
0093                 //    (make sure to handle rollover of phi properly)
0094                 const double de_rel = std::abs((pc.getEnergy() - ec.getEnergy()) / ec.getEnergy());
0095                 const double deta = std::abs(edm4hep::utils::eta(pc.getPosition())
0096                                            - edm4hep::utils::eta(ec.getPosition()));
0097                 // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow
0098                 // for phi rollovers
0099                 const double dphi = edm4hep::utils::angleAzimuthal(pc.getPosition())
0100                                   - edm4hep::utils::angleAzimuthal(ec.getPosition());
0101                 const double dsphi = std::abs(sin(0.5 * dphi));
0102                 if ((m_cfg.energyRelTolerance > 0 && de_rel > m_cfg.energyRelTolerance) ||
0103                     (m_cfg.etaTolerance > 0 && deta > m_cfg.etaTolerance) ||
0104                     (m_cfg.phiTolerance > 0 && dsphi > sin(0.5 * m_cfg.phiTolerance))) {
0105                     continue;
0106                 }
0107                 // --> if we get here, we have a match within tolerance. Now treat the case
0108                 //     where we have multiple matches. In this case take the one with the closest
0109                 //     energies.
0110                 // 2. best match?
0111                 const double delta = fabs(pc.getEnergy() - ec.getEnergy());
0112                 if (delta < best_delta) {
0113                     best_delta = delta;
0114                     best_match = ie;
0115                 }
0116             }
0117
0118             // Create a merged cluster if we find a good match
0119             if (best_match >= 0) {
0120
0121                 const auto& ec = (*energy_clus)[best_match];
0122
0123                 auto new_clus  = merged_clus->create();
0124                 new_clus.setEnergy(ec.getEnergy());
0125                 new_clus.setEnergyError(ec.getEnergyError());
0126                 new_clus.setTime(pc.getTime());
0127                 new_clus.setNhits(pc.getNhits() + ec.getNhits());
0128                 new_clus.setPosition(pc.getPosition());
0129                 new_clus.setPositionError(pc.getPositionError());
0130                 new_clus.addToClusters(pc);
0131                 new_clus.addToClusters(ec);
0132
0133                 trace("   --> Found matching energy cluster {}, energy: {}", ec.getObjectID().index, ec.getEnergy() );
0134                 trace("   --> Created a new combined cluster {}, energy: {}", new_clus.getObjectID().index, new_clus.getEnergy() );
0135
0136                 // find association from energy cluster
0137                 auto ea = energy_assoc->begin();
0138                 for (; ea != energy_assoc->end(); ++ea) {
0139                     if (ea->getRec() == ec) {
0140                         break;
0141                     }
0142                 }
0143                 // find association from position cluster if different
0144                 auto pa = pos_assoc->begin();
0145                 for (; pa != pos_assoc->end(); ++pa) {
0146                     if (pa->getRec() == pc) {
0147                         break;
0148                     }
0149                 }
0150                 if (ea != energy_assoc->end() || pa != pos_assoc->end()) {
0151                     // we must write an association
0152                     if (ea != energy_assoc->end() && pa != pos_assoc->end()) {
0153                         // we have two associations
0154                         if (pa->getSimID() == ea->getSimID()) {
0155                             // both associations agree on the MCParticles entry
0156                             auto clusterassoc = merged_assoc->create();
0157                             clusterassoc.setRecID(new_clus.getObjectID().index);
0158                             clusterassoc.setSimID(ea->getSimID());
0159                             clusterassoc.setWeight(1.0);
0160                             clusterassoc.setRec(new_clus);
0161                             clusterassoc.setSim(ea->getSim());
0162                         } else {
0163                             // both associations disagree on the MCParticles entry
0164                             debug("   --> Two associations added to {} and {}", ea->getSimID(), pa->getSimID());
0165                             auto clusterassoc1 = merged_assoc->create();
0166                             clusterassoc1.setRecID(new_clus.getObjectID().index);
0167                             clusterassoc1.setSimID(ea->getSimID());
0168                             clusterassoc1.setWeight(0.5);
0169                             clusterassoc1.setRec(new_clus);
0170                             clusterassoc1.setSim(ea->getSim());
0171                             auto clusterassoc2 = merged_assoc->create();
0172                             clusterassoc2.setRecID(new_clus.getObjectID().index);
0173                             clusterassoc2.setSimID(pa->getSimID());
0174                             clusterassoc2.setWeight(0.5);
0175                             clusterassoc2.setRec(new_clus);
0176                             clusterassoc2.setSim(pa->getSim());
0177                         }
0178                     } else if (ea != energy_assoc->end()) {
0179                         // no position association
0180                         debug("   --> Only added energy cluster association to {}", ea->getSimID());
0181                         auto clusterassoc = merged_assoc->create();
0182                         clusterassoc.setRecID(new_clus.getObjectID().index);
0183                         clusterassoc.setSimID(ea->getSimID());
0184                         clusterassoc.setWeight(1.0);
0185                         clusterassoc.setRec(new_clus);
0186                         clusterassoc.setSim(ea->getSim());
0187                     } else if (pa != pos_assoc->end()) {
0188                         // no energy association
0189                         debug("   --> Only added position cluster association to {}", pa->getSimID());
0190                         auto clusterassoc = merged_assoc->create();
0191                         clusterassoc.setRecID(new_clus.getObjectID().index);
0192                         clusterassoc.setSimID(pa->getSimID());
0193                         clusterassoc.setWeight(1.0);
0194                         clusterassoc.setRec(new_clus);
0195                         clusterassoc.setSim(pa->getSim());
0196                     }
0197                 }
0198
0199                 // label our energy cluster as consumed
0200                 consumed[best_match] = true;
0201
0202                 debug("  Matched position cluster {} with energy cluster {}", pc.getObjectID().index, ec.getObjectID().index);
0203                 debug("  - Position cluster: (E: {}, phi: {}, z: {})", pc.getEnergy(),
0204                              edm4hep::utils::angleAzimuthal(pc.getPosition()), pc.getPosition().z);
0205                 debug("  - Energy cluster: (E: {}, phi: {}, z: {})", ec.getEnergy(),
0206                              edm4hep::utils::angleAzimuthal(ec.getPosition()), ec.getPosition().z);
0207                 debug("  ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.getEnergy(),
0208                              edm4hep::utils::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z);
0209
0210
0211             } else {
0212
0213                 debug("  Unmatched position cluster {}", pc.getObjectID().index);
0214
0215             }
0216
0217         }
0218     }
0219   };
0220
0221 } // namespace eicrecon