algorithms/calorimetry/CalorimeterTruthClustering.cc

0001 // Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
0002 // SPDX-License-Identifier: LGPL-3.0-or-later
0003
0004 #include "CalorimeterTruthClustering.h"
0005
0006 #include <DD4hep/config.h>
0007 #include <edm4hep/CaloHitContributionCollection.h>
0008 #include <edm4hep/MCParticleCollection.h>
0009 #include <podio/ObjectID.h>
0010 #include <cstddef>
0011 #include <cstdint>
0012 #include <gsl/pointers>
0013 #include <map>
0014
0015 using namespace dd4hep;
0016
0017 namespace eicrecon {
0018
0019   void CalorimeterTruthClustering::init() { }
0020
0021
0022   void CalorimeterTruthClustering::process(
0023       const CalorimeterTruthClustering::Input& input,
0024       const CalorimeterTruthClustering::Output& output) const {
0025     const auto [hits, mc] = input;
0026     auto [clusters] = output;
0027
0028     // Map mc track ID to protoCluster index
0029     std::map<int32_t, int32_t> protoIndex;
0030
0031     // Loop over all calorimeter hits and sort per mcparticle
0032     for (const auto& hit : *hits) {
0033         // The original algorithm used the following to get the mcHit:
0034         //
0035         //    const auto& mcHit     = mc[hit->getObjectID().index];
0036         //
0037         // This assumes there is a one-to-one relation between the truth hit
0038         // (hits) and the reconstructed hit (mc). At least insofar as the
0039         // index in "hits" is being used to index the "mc" container.
0040         //
0041         // If the objects in "hits" have not been added to a collection,
0042         // then they will have getObjectID().index = "untracked" = -1
0043         //
0044         // The way we handle this is here is to check if getObjectID().index
0045         // is within the size limits of mc which includes >=0. If so, then
0046         // assume the old code is valid. If not, then we need to search
0047         // for the right hit.
0048         // FIXME: This is clearly not the right way to do this! Podio needs
0049         // FIXME: to be fixed so proper object tracking can be done without
0050         // FIXME: requiring Collection classes be used to manage all objects.
0051         std::size_t mcIndex;
0052         if ((hit.getObjectID().index >= 0) && (hit.getObjectID().index < mc->size())) {
0053             mcIndex = hit.getObjectID().index;
0054         } else {
0055             mcIndex = 0;
0056             bool success = false;
0057             for (auto tmpmc : *mc) {
0058                 if (tmpmc.getCellID() == hit.getCellID()) {
0059                     success = true;
0060                     break;
0061                 }
0062                 mcIndex++;
0063             }
0064             if (not success) {
0065                 continue; // ignore hit if we couldn't match it to truth hit
0066             }
0067         }
0068
0069         const auto &trackID = (*mc)[mcIndex].getContributions(0).getParticle().getObjectID().index;
0070         // Create a new protocluster if we don't have one for this trackID
0071         if (protoIndex.count(trackID) == 0) {
0072             clusters->create();
0073             protoIndex[trackID] = clusters->size() - 1;
0074         }
0075         // Add hit to the appropriate protocluster
0076         (*clusters)[protoIndex[trackID]].addToHits(hit);
0077         (*clusters)[protoIndex[trackID]].addToWeights(1);
0078     }
0079
0080   }
0081
0082 } // namespace eicrecon