calorimetry/todo/CalorimeterHitsMerger.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
0003
0004 /*
0005  *  An algorithm to group readout hits from a calorimeter
0006  *  Energy is summed
0007  *
0008  *  Author: Chao Peng (ANL), 03/31/2021
0009  */
0010 #include <algorithm>
0011 #include <bitset>
0012 #include <tuple>
0013 #include <unordered_map>
0014
0015 #include "Gaudi/Property.h"
0016 #include "GaudiAlg/GaudiAlgorithm.h"
0017 #include "GaudiAlg/GaudiTool.h"
0018 #include "GaudiAlg/Transformer.h"
0019 #include "GaudiKernel/PhysicalConstants.h"
0020 #include "GaudiKernel/RndmGenerators.h"
0021 #include "GaudiKernel/ToolHandle.h"
0022
0023 #include "DDRec/CellIDPositionConverter.h"
0024 #include "DDRec/Surface.h"
0025 #include "DDRec/SurfaceManager.h"
0026 #include "DDSegmentation/BitFieldCoder.h"
0027
0028 #include "fmt/format.h"
0029 #include "fmt/ranges.h"
0030
0031 #include "JugBase/DataHandle.h"
0032 #include "JugBase/IGeoSvc.h"
0033
0034 // Event Model related classes
0035 #include "edm4eic/CalorimeterHitCollection.h"
0036
0037 using namespace Gaudi::Units;
0038
0039 namespace Jug::Reco {
0040
0041 /** Calorimeter hits merging algorithm.
0042  *
0043  *  An algorithm to group readout hits from a calorimeter
0044  *  Energy is summed
0045  *
0046  *  \ingroup reco
0047  */
0048 class CalorimeterHitsMerger : public GaudiAlgorithm {
0049 private:
0050   Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0051   Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
0052   // field names to generate id mask, the hits will be grouped by masking the field
0053   Gaudi::Property<std::vector<std::string>> u_fields{this, "fields", {"layer"}};
0054   // reference field numbers to locate position for each merged hits group
0055   Gaudi::Property<std::vector<int>> u_refs{this, "fieldRefNumbers", {}};
0056   DataHandle<edm4eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0057   DataHandle<edm4eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0058                                                                   this};
0059
0060   SmartIF<IGeoSvc> m_geoSvc;
0061   uint64_t id_mask{0}, ref_mask{0};
0062
0063 public:
0064   CalorimeterHitsMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0065     declareProperty("inputHitCollection", m_inputHitCollection, "");
0066     declareProperty("outputHitCollection", m_outputHitCollection, "");
0067   }
0068
0069   StatusCode initialize() override {
0070     if (GaudiAlgorithm::initialize().isFailure()) {
0071       return StatusCode::FAILURE;
0072     }
0073
0074     m_geoSvc = service(m_geoSvcName);
0075     if (!m_geoSvc) {
0076       error() << "Unable to locate Geometry Service. "
0077               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0078       return StatusCode::FAILURE;
0079     }
0080
0081     if (m_readout.value().empty()) {
0082       error() << "readoutClass is not provided, it is needed to know the fields in readout ids" << endmsg;
0083       return StatusCode::FAILURE;
0084     }
0085
0086     try {
0087       auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
0088       id_mask      = 0;
0089       std::vector<std::pair<std::string, int>> ref_fields;
0090       for (size_t i = 0; i < u_fields.size(); ++i) {
0091         id_mask |= id_desc.field(u_fields[i])->mask();
0092         // use the provided id number to find ref cell, or use 0
0093         int ref = i < u_refs.size() ? u_refs[i] : 0;
0094         ref_fields.emplace_back(u_fields[i], ref);
0095       }
0096       ref_mask = id_desc.encode(ref_fields);
0097       // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
0098     } catch (...) {
0099       error() << "Failed to load ID decoder for " << m_readout << endmsg;
0100       return StatusCode::FAILURE;
0101     }
0102     id_mask = ~id_mask;
0103     info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg;
0104     return StatusCode::SUCCESS;
0105   }
0106
0107   StatusCode execute() override {
0108     // input collections
0109     const auto& inputs = *m_inputHitCollection.get();
0110     // Create output collections
0111     auto& outputs = *m_outputHitCollection.createAndPut();
0112
0113     // find the hits that belong to the same group (for merging)
0114     std::unordered_map<long long, std::vector<edm4eic::CalorimeterHit>> merge_map;
0115     for (const auto& h : inputs) {
0116       int64_t id = h.getCellID() & id_mask;
0117       // use the reference field position
0118       auto it = merge_map.find(id);
0119       if (it == merge_map.end()) {
0120         merge_map[id] = {h};
0121       } else {
0122         it->second.push_back(h);
0123       }
0124     }
0125
0126     // sort hits by energy from large to small
0127     std::for_each(merge_map.begin(), merge_map.end(), [](auto& it) {
0128       std::sort(it.second.begin(), it.second.end(), [](const auto& h1, const auto& h2) {
0129         return h1.getEnergy() > h2.getEnergy();
0130       });
0131     });
0132
0133     // reconstruct info for merged hits
0134     // dd4hep decoders
0135     auto poscon = m_geoSvc->cellIDPositionConverter();
0136     auto volman = m_geoSvc->detector()->volumeManager();
0137
0138     for (auto& [id, hits] : merge_map) {
0139       // reference fields id
0140       const uint64_t ref_id = id | ref_mask;
0141       // global positions
0142       const auto gpos = poscon->position(ref_id);
0143       // local positions
0144       auto alignment = volman.lookupDetElement(ref_id).nominal();
0145       const auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0146       debug() << volman.lookupDetElement(ref_id).path() << ", " << volman.lookupDetector(ref_id).path() << endmsg;
0147       // sum energy
0148       float energy      = 0.;
0149       float energyError = 0.;
0150       float time        = 0;
0151       float timeError   = 0;
0152       for (auto& hit : hits) {
0153         energy += hit.getEnergy();
0154         energyError += hit.getEnergyError() * hit.getEnergyError();
0155         time += hit.getTime();
0156         timeError += hit.getTimeError() * hit.getTimeError();
0157       }
0158       energyError = sqrt(energyError);
0159       time /= hits.size();
0160       timeError = sqrt(timeError) / hits.size();
0161
0162       const auto& href = hits.front();
0163
0164       // create const vectors for passing to hit initializer list
0165       const decltype(edm4eic::CalorimeterHitData::position) position(
0166         gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm
0167       );
0168       const decltype(edm4eic::CalorimeterHitData::local) local(
0169         pos.x(), pos.y(), pos.z()
0170       );
0171
0172       outputs.push_back(
0173           edm4eic::CalorimeterHit{href.getCellID(),
0174                               energy,
0175                               energyError,
0176                               time,
0177                               timeError,
0178                               position,
0179                               href.getDimension(),
0180                               href.getSector(),
0181                               href.getLayer(),
0182                               local}); // Can do better here? Right now position is mapped on the central hit
0183     }
0184
0185     if (msgLevel(MSG::DEBUG)) {
0186       debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
0187     }
0188
0189     return StatusCode::SUCCESS;
0190   }
0191
0192 }; // class CalorimeterHitsMerger
0193
0194 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0195 DECLARE_COMPONENT(CalorimeterHitsMerger)
0196
0197 } // namespace Jug::Reco