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 void CalorimeterTruthClustering::process(const CalorimeterTruthClustering::Input& input,
0022                                          const CalorimeterTruthClustering::Output& output) const {
0023   const auto [hits, mc] = input;
0024   auto [clusters]       = output;
0025
0026   // Map mc track ID to protoCluster index
0027   std::map<int32_t, int32_t> protoIndex;
0028
0029   // Loop over all calorimeter hits and sort per mcparticle
0030   for (const auto& hit : *hits) {
0031     // The original algorithm used the following to get the mcHit:
0032     //
0033     //    const auto& mcHit     = mc[hit->getObjectID().index];
0034     //
0035     // This assumes there is a one-to-one relation between the truth hit
0036     // (hits) and the reconstructed hit (mc). At least insofar as the
0037     // index in "hits" is being used to index the "mc" container.
0038     //
0039     // If the objects in "hits" have not been added to a collection,
0040     // then they will have getObjectID().index = "untracked" = -1
0041     //
0042     // The way we handle this is here is to check if getObjectID().index
0043     // is within the size limits of mc which includes >=0. If so, then
0044     // assume the old code is valid. If not, then we need to search
0045     // for the right hit.
0046     // FIXME: This is clearly not the right way to do this! Podio needs
0047     // FIXME: to be fixed so proper object tracking can be done without
0048     // FIXME: requiring Collection classes be used to manage all objects.
0049     std::size_t mcIndex = 0;
0050     if ((hit.getObjectID().index >= 0) &&
0051         (hit.getObjectID().index < static_cast<long>(mc->size()))) {
0052       mcIndex = hit.getObjectID().index;
0053     } else {
0054       mcIndex      = 0;
0055       bool success = false;
0056       for (auto tmpmc : *mc) {
0057         if (tmpmc.getCellID() == hit.getCellID()) {
0058           success = true;
0059           break;
0060         }
0061         mcIndex++;
0062       }
0063       if (not success) {
0064         continue; // ignore hit if we couldn't match it to truth hit
0065       }
0066     }
0067
0068     const auto& trackID = (*mc)[mcIndex].getContributions(0).getParticle().getObjectID().index;
0069     // Create a new protocluster if we don't have one for this trackID
0070     if (!protoIndex.contains(trackID)) {
0071       clusters->create();
0072       protoIndex[trackID] = clusters->size() - 1;
0073     }
0074     // Add hit to the appropriate protocluster
0075     (*clusters)[protoIndex[trackID]].addToHits(hit);
0076     (*clusters)[protoIndex[trackID]].addToWeights(1);
0077   }
0078 }
0079
0080 } // namespace eicrecon