File indexing completed on 2025-08-06 08:10:48
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "ActsExamples/Io/EDM4hep/EDM4hepReader.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0015 #include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp"
0016 #include "ActsExamples/EventData/SimHit.hpp"
0017 #include "ActsExamples/EventData/SimParticle.hpp"
0018 #include "ActsExamples/Framework/WhiteBoard.hpp"
0019 #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp"
0020 #include "ActsExamples/Utilities/Paths.hpp"
0021 #include "ActsFatras/EventData/Barcode.hpp"
0022
0023 #include <algorithm>
0024 #include <iomanip>
0025 #include <map>
0026 #include <stdexcept>
0027
0028 #include <edm4hep/MCParticle.h>
0029 #include <edm4hep/SimTrackerHit.h>
0030 #include <edm4hep/SimTrackerHitCollection.h>
0031 #include <podio/Frame.h>
0032 #include <podio/ObjectID.h>
0033 #include <podio/ROOTFrameReader.h>
0034
0035 namespace ActsExamples {
0036
0037 EDM4hepReader::EDM4hepReader(const Config& config, Acts::Logging::Level level)
0038 : m_cfg(config),
0039 m_logger(Acts::getDefaultLogger("EDM4hepParticleReader", level)) {
0040 if (m_cfg.outputParticlesInitial.empty()) {
0041 throw std::invalid_argument("Missing output collection initial particles");
0042 }
0043
0044 if (m_cfg.outputParticlesFinal.empty()) {
0045 throw std::invalid_argument("Missing output collection final particles");
0046 }
0047
0048 if (m_cfg.outputParticlesGenerator.empty()) {
0049 throw std::invalid_argument(
0050 "Missing output collection generator particles");
0051 }
0052
0053 if (m_cfg.outputSimHits.empty()) {
0054 throw std::invalid_argument("Missing output collection sim hits");
0055 }
0056
0057 m_eventsRange = std::make_pair(0, reader().getEntries("events"));
0058
0059 m_outputParticlesInitial.initialize(m_cfg.outputParticlesInitial);
0060 m_outputParticlesFinal.initialize(m_cfg.outputParticlesFinal);
0061 m_outputParticlesGenerator.initialize(m_cfg.outputParticlesGenerator);
0062 m_outputSimHits.initialize(m_cfg.outputSimHits);
0063
0064 m_cfg.trackingGeometry->visitSurfaces([&](const auto* surface) {
0065 const auto* detElement = dynamic_cast<const Acts::DD4hepDetectorElement*>(
0066 surface->associatedDetectorElement());
0067
0068 if (detElement == nullptr) {
0069 ACTS_ERROR("Surface has no associated detector element");
0070 return;
0071 }
0072
0073 const auto translation = detElement->sourceElement()
0074 .nominal()
0075 .worldTransformation()
0076 .GetTranslation();
0077 Acts::Vector3 position;
0078 position << translation[0], translation[1], translation[2];
0079 position *= Acts::UnitConstants::cm;
0080
0081 m_surfaceMap.insert({detElement->sourceElement().key(), surface});
0082 });
0083 }
0084
0085 podio::ROOTFrameReader& EDM4hepReader::reader() {
0086 bool exists = false;
0087 auto& reader = m_reader.local(exists);
0088 if (!exists) {
0089 reader.openFile(m_cfg.inputPath);
0090 }
0091
0092 return reader;
0093 }
0094
0095 std::string EDM4hepReader::name() const {
0096 return "EDM4hepReader";
0097 }
0098
0099 std::pair<std::size_t, std::size_t> EDM4hepReader::availableEvents() const {
0100 return m_eventsRange;
0101 }
0102
0103 namespace {
0104 std::string vid(unsigned int vtx) {
0105 return "V" + std::to_string(vtx);
0106 }
0107
0108 std::string pid(const SimParticle& particle) {
0109 return "P" + std::to_string(particle.particleId().value());
0110 }
0111
0112 std::string plabel(const SimParticle& particle) {
0113 using namespace Acts::UnitLiterals;
0114 std::stringstream ss;
0115 ss << particle.pdg() << "\\n(" << particle.particleId() << ")\\n"
0116 << "p=" << std::setprecision(3) << particle.absoluteMomentum() / 1_GeV
0117 << " GeV";
0118 return ss.str();
0119 }
0120
0121 }
0122
0123 void EDM4hepReader::graphviz(
0124 std::ostream& os, const SimParticleContainer::sequence_type& particles,
0125 const ParentRelationship& parents) const {
0126 os << "digraph Event {\n";
0127
0128 std::set<unsigned int> primaryVertices;
0129
0130 for (const auto& particle : particles) {
0131 if (particle.particleId().generation() == 0) {
0132 primaryVertices.insert(particle.particleId().vertexPrimary());
0133
0134 os << vid(particle.particleId().vertexPrimary()) << " -> "
0135 << pid(particle) << ";\n";
0136 }
0137
0138 os << pid(particle) << " [label=\"" << plabel(particle) << "\"];\n";
0139 }
0140
0141 for (const auto [childIdx, parentIdx] : parents) {
0142 const auto& child = particles[childIdx];
0143 const auto& parent = particles[parentIdx];
0144 os << pid(parent) << " -> " << pid(child);
0145
0146 if (parent.particleId().vertexSecondary() ==
0147 child.particleId().vertexSecondary()) {
0148 os << " [style=dashed]";
0149 }
0150
0151 os << ";\n";
0152 }
0153
0154 for (unsigned int vtx : primaryVertices) {
0155 os << vid(vtx) << " [label=\"PV" << vtx << "\"];\n";
0156 }
0157
0158 os << "}";
0159 }
0160
0161 ProcessCode EDM4hepReader::read(const AlgorithmContext& ctx) {
0162 podio::Frame frame = reader().readEntry("events", ctx.eventNumber);
0163 const auto& mcParticleCollection =
0164 frame.get<edm4hep::MCParticleCollection>(m_cfg.inputParticles);
0165
0166 ACTS_DEBUG("Reading EDM4hep inputs");
0167
0168 SimParticleContainer::sequence_type unordered;
0169
0170
0171
0172
0173 std::vector<std::pair<Acts::Vector3, std::vector<edm4hep::MCParticle>>>
0174 primaryVertices;
0175 for (const auto& particle : mcParticleCollection) {
0176 if (particle.parents_size() > 0) {
0177
0178 continue;
0179 }
0180 const auto& vtx = particle.getVertex();
0181 Acts::Vector3 vtxPos = {vtx[0], vtx[1], vtx[2]};
0182 vtxPos /= Acts::UnitConstants::mm;
0183
0184
0185 auto it = std::find_if(
0186 primaryVertices.begin(), primaryVertices.end(),
0187 [&vtxPos](
0188 const std::pair<Acts::Vector3, std::vector<edm4hep::MCParticle>>&
0189 pair) { return pair.first == vtxPos; });
0190
0191 if (it == primaryVertices.end()) {
0192 ACTS_DEBUG("Found primary vertex at " << vtx.x << ", " << vtx.y << ", "
0193 << vtx.z);
0194 primaryVertices.push_back({vtxPos, {particle}});
0195 } else {
0196 it->second.push_back(particle);
0197 }
0198 }
0199
0200 ACTS_DEBUG("Found " << primaryVertices.size() << " primary vertices");
0201
0202
0203 ParentRelationship parentRelationship;
0204
0205
0206
0207 std::unordered_map<int, std::size_t> edm4hepParticleMap;
0208
0209 std::size_t nPrimaryVertices = 0;
0210
0211 for (const auto& [vtxPos, particles] : primaryVertices) {
0212 nPrimaryVertices += 1;
0213 ACTS_DEBUG("Walking particle tree for primary vertex at "
0214 << vtxPos.x() << ", " << vtxPos.y() << ", " << vtxPos.z());
0215 std::size_t nParticles = 0;
0216 std::size_t nSecondaryVertices = 0;
0217 std::size_t maxGen = 0;
0218 auto startSize = unordered.size();
0219 for (const auto& inParticle : particles) {
0220 nParticles += 1;
0221 SimParticle particle{EDM4hepUtil::readParticle(inParticle)};
0222 particle.setParticleId(SimBarcode{}
0223 .setParticle(nParticles)
0224 .setVertexPrimary(nPrimaryVertices));
0225 ACTS_VERBOSE("+ add particle " << particle);
0226 ACTS_VERBOSE(" - at " << particle.position().transpose());
0227 ACTS_VERBOSE(" - createdInSim: " << inParticle.isCreatedInSimulation());
0228 ACTS_VERBOSE(" - vertexIsNotEndpointOfParent: "
0229 << inParticle.vertexIsNotEndpointOfParent());
0230 ACTS_VERBOSE(" - isStopped: " << inParticle.isStopped());
0231 ACTS_VERBOSE(" - endpoint: " << inParticle.getEndpoint().x << ", "
0232 << inParticle.getEndpoint().y << ", "
0233 << inParticle.getEndpoint().z);
0234 const auto pid = particle.particleId();
0235 unordered.push_back(std::move(particle));
0236 edm4hepParticleMap[inParticle.getObjectID().index] = unordered.size() - 1;
0237 processChildren(inParticle, pid, unordered, parentRelationship,
0238 edm4hepParticleMap, nSecondaryVertices, maxGen);
0239 }
0240 ACTS_VERBOSE("Primary vertex complete, produced "
0241 << (unordered.size() - startSize) << " particles and "
0242 << nSecondaryVertices << " secondary vertices in " << maxGen
0243 << " generations");
0244 setSubParticleIds(std::next(unordered.begin(), startSize), unordered.end());
0245 }
0246
0247 ACTS_DEBUG("Found " << unordered.size() << " particles");
0248
0249
0250
0251 SimParticleContainer particlesFinal;
0252 SimParticleContainer particlesGenerator;
0253 for (const auto& inParticle : mcParticleCollection) {
0254 const std::size_t index =
0255 edm4hepParticleMap.find(inParticle.getObjectID().index)->second;
0256 const auto& particleInitial = unordered.at(index);
0257 if (!inParticle.isCreatedInSimulation()) {
0258 particlesGenerator.insert(particleInitial);
0259 }
0260 SimParticle particleFinal = particleInitial;
0261
0262 float time = inParticle.getTime() * Acts::UnitConstants::ns;
0263 for (const auto& daughter : inParticle.getDaughters()) {
0264 if (!daughter.vertexIsNotEndpointOfParent()) {
0265 time = daughter.getTime() * Acts::UnitConstants::ns;
0266 break;
0267 }
0268 }
0269
0270 particleFinal.setPosition4(
0271 inParticle.getEndpoint()[0] * Acts::UnitConstants::mm,
0272 inParticle.getEndpoint()[1] * Acts::UnitConstants::mm,
0273 inParticle.getEndpoint()[2] * Acts::UnitConstants::mm, time);
0274
0275 Acts::Vector3 momentumFinal = {inParticle.getMomentumAtEndpoint()[0],
0276 inParticle.getMomentumAtEndpoint()[1],
0277 inParticle.getMomentumAtEndpoint()[2]};
0278 particleFinal.setDirection(momentumFinal.normalized());
0279 particleFinal.setAbsoluteMomentum(momentumFinal.norm());
0280
0281 ACTS_VERBOSE("- Updated particle initial -> final, position: "
0282 << particleInitial.fourPosition().transpose() << " -> "
0283 << particleFinal.fourPosition().transpose());
0284 ACTS_VERBOSE(" momentum: "
0285 << particleInitial.fourMomentum().transpose() << " -> "
0286 << particleFinal.fourMomentum().transpose());
0287
0288 particlesFinal.insert(particleFinal);
0289 }
0290
0291
0292 SimParticleContainer particlesInitial;
0293 particlesInitial.insert(unordered.begin(), unordered.end());
0294
0295 if (!m_cfg.graphvizOutput.empty()) {
0296 std::string path = perEventFilepath(m_cfg.graphvizOutput, "particles.dot",
0297 ctx.eventNumber);
0298 std::ofstream dot(path);
0299 graphviz(dot, unordered, parentRelationship);
0300 }
0301
0302 SimHitContainer simHits;
0303
0304 ACTS_DEBUG("Reading sim hits from " << m_cfg.inputSimHits.size()
0305 << " sim hit collections");
0306 for (const auto& name : m_cfg.inputSimHits) {
0307 const auto& inputHits = frame.get<edm4hep::SimTrackerHitCollection>(name);
0308
0309 for (const auto& hit : inputHits) {
0310 auto simHit = EDM4hepUtil::readSimHit(
0311 hit,
0312 [&](const auto& inParticle) {
0313 ACTS_VERBOSE("SimHit has source particle: "
0314 << hit.getMCParticle().getObjectID().index);
0315 auto it = edm4hepParticleMap.find(inParticle.getObjectID().index);
0316 if (it == edm4hepParticleMap.end()) {
0317 ACTS_ERROR(
0318 "SimHit has source particle that we did not see before");
0319 return SimBarcode{};
0320 }
0321 const auto& particle = unordered.at(it->second);
0322 ACTS_VERBOSE("- " << inParticle.getObjectID().index << " -> "
0323 << particle.particleId());
0324 return particle.particleId();
0325 },
0326 [&](std::uint64_t cellId) {
0327 ACTS_VERBOSE("CellID: " << cellId);
0328
0329 const auto& vm = m_cfg.dd4hepDetector->geometryService->detector()
0330 .volumeManager();
0331
0332 const auto detElement = vm.lookupDetElement(cellId);
0333
0334 ACTS_VERBOSE(" -> detElement: " << detElement.name());
0335 ACTS_VERBOSE(" -> id: " << detElement.id());
0336 ACTS_VERBOSE(" -> key: " << detElement.key());
0337
0338 Acts::Vector3 position;
0339 position << detElement.nominal()
0340 .worldTransformation()
0341 .GetTranslation()[0],
0342 detElement.nominal().worldTransformation().GetTranslation()[1],
0343 detElement.nominal().worldTransformation().GetTranslation()[2];
0344 position *= Acts::UnitConstants::cm;
0345
0346 ACTS_VERBOSE(" -> detElement position: " << position.transpose());
0347
0348 auto it = m_surfaceMap.find(detElement.key());
0349 if (it == m_surfaceMap.end()) {
0350 ACTS_ERROR("Unable to find surface for detElement "
0351 << detElement.name() << " with cellId " << cellId);
0352 }
0353 const auto* surface = it->second;
0354 ACTS_VERBOSE(" -> surface: " << surface->geometryId());
0355 return surface->geometryId();
0356 });
0357
0358 simHits.insert(std::move(simHit));
0359 }
0360 }
0361
0362 if (m_cfg.sortSimHitsInTime) {
0363 ACTS_DEBUG("Sorting sim hits in time");
0364 std::multimap<ActsFatras::Barcode, std::size_t> hitsByParticle;
0365
0366 for (std::size_t i = 0; i < simHits.size(); ++i) {
0367 hitsByParticle.insert({simHits.nth(i)->particleId(), i});
0368 }
0369
0370 for (auto it = hitsByParticle.begin(), end = hitsByParticle.end();
0371 it != end; it = hitsByParticle.upper_bound(it->first)) {
0372 std::cout << "Particle " << it->first << " has "
0373 << hitsByParticle.count(it->first) << " hits" << std::endl;
0374
0375 std::vector<std::size_t> hitIndices;
0376 hitIndices.reserve(hitsByParticle.count(it->first));
0377 for (auto hitIndex : makeRange(hitsByParticle.equal_range(it->first))) {
0378 hitIndices.push_back(hitIndex.second);
0379 }
0380
0381 if (logger().doPrint(Acts::Logging::VERBOSE)) {
0382 ACTS_VERBOSE("Before sorting:");
0383 for (const auto& hitIdx : hitIndices) {
0384 ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index()
0385 << " " << simHits.nth(hitIdx)->time());
0386 }
0387 }
0388
0389 std::sort(hitIndices.begin(), hitIndices.end(),
0390 [&](std::size_t a, std::size_t b) {
0391 return simHits.nth(a)->time() < simHits.nth(b)->time();
0392 });
0393
0394 for (std::size_t i = 0; i < hitIndices.size(); ++i) {
0395 auto& hit = *simHits.nth(hitIndices[i]);
0396 SimHit updatedHit{hit.geometryId(), hit.particleId(),
0397 hit.fourPosition(), hit.momentum4Before(),
0398 hit.momentum4After(), int32_t(i)};
0399 hit = updatedHit;
0400 }
0401
0402 if (logger().doPrint(Acts::Logging::VERBOSE)) {
0403 ACTS_VERBOSE("After sorting:");
0404 for (const auto& hitIdx : hitIndices) {
0405 ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index()
0406 << " " << simHits.nth(hitIdx)->time());
0407 }
0408 }
0409 }
0410 }
0411
0412 m_outputParticlesInitial(ctx, std::move(particlesInitial));
0413 m_outputParticlesFinal(ctx, std::move(particlesFinal));
0414 m_outputParticlesGenerator(ctx, std::move(particlesGenerator));
0415
0416 m_outputSimHits(ctx, std::move(simHits));
0417
0418 return ProcessCode::SUCCESS;
0419 }
0420
0421 void EDM4hepReader::processChildren(
0422 const edm4hep::MCParticle& inParticle, SimBarcode parentId,
0423 SimParticleContainer::sequence_type& particles,
0424 ParentRelationship& parentRelationship,
0425 std::unordered_map<int, std::size_t>& particleMap,
0426 std::size_t& nSecondaryVertices, std::size_t& maxGen) const {
0427 constexpr auto indent = [&](std::size_t n) {
0428 std::string result;
0429 for (std::size_t i = 0; i < n; ++i) {
0430 result += " ";
0431 }
0432 return result;
0433 };
0434
0435 const std::size_t gen = parentId.generation();
0436 maxGen = std::max(maxGen, gen);
0437
0438 ACTS_VERBOSE(indent(gen) << " - processing daughters for input particle "
0439 << inParticle.id());
0440 ACTS_VERBOSE(indent(gen) << " -> found " << inParticle.daughters_size()
0441 << " daughter(s)");
0442
0443 bool parentDecayed =
0444 std::any_of(inParticle.daughters_begin(), inParticle.daughters_end(),
0445 [](const edm4hep::MCParticle& daughter) {
0446 return !daughter.vertexIsNotEndpointOfParent();
0447 });
0448 std::size_t secondaryVertex = 0;
0449 if (parentDecayed) {
0450 ACTS_VERBOSE(indent(gen) << " -> parent decays");
0451 secondaryVertex = ++nSecondaryVertices;
0452 }
0453
0454 std::size_t parentIndex = particles.size() - 1;
0455
0456 std::size_t nParticles = 0;
0457 for (const auto& daughter : inParticle.getDaughters()) {
0458 SimParticle particle = EDM4hepUtil::readParticle(daughter);
0459
0460 auto pid = parentId.makeDescendant(nParticles++);
0461 if (daughter.vertexIsNotEndpointOfParent()) {
0462
0463 } else {
0464
0465 pid = pid.setVertexSecondary(secondaryVertex);
0466 }
0467 particle.setParticleId(pid);
0468
0469 ACTS_VERBOSE(indent(particle.particleId().generation())
0470 << "+ add particle " << particle);
0471 ACTS_VERBOSE(indent(particle.particleId().generation())
0472 << " - generation: " << particle.particleId().generation());
0473 ACTS_VERBOSE(indent(particle.particleId().generation())
0474 << " - at " << particle.position().transpose());
0475 ACTS_VERBOSE(indent(particle.particleId().generation())
0476 << " - createdInSim: "
0477 << daughter.isCreatedInSimulation());
0478 ACTS_VERBOSE(indent(particle.particleId().generation())
0479 << " - vertexIsNotEndpointOfParent: "
0480 << daughter.vertexIsNotEndpointOfParent());
0481 ACTS_VERBOSE(indent(particle.particleId().generation())
0482 << " - isStopped: " << daughter.isStopped());
0483 ACTS_VERBOSE(indent(particle.particleId().generation())
0484 << " - endpoint: " << daughter.getEndpoint().x << ", "
0485 << daughter.getEndpoint().y << ", "
0486 << daughter.getEndpoint().z);
0487
0488 particles.push_back(std::move(particle));
0489 particleMap[daughter.getObjectID().index] = particles.size() - 1;
0490 parentRelationship[particles.size() - 1] = parentIndex;
0491 processChildren(daughter, pid, particles, parentRelationship, particleMap,
0492 nSecondaryVertices, maxGen);
0493 }
0494 }
0495
0496 void EDM4hepReader::setSubParticleIds(
0497 const SimParticleContainer::sequence_type::iterator& begin,
0498 const SimParticleContainer::sequence_type::iterator& end) {
0499 std::vector<std::size_t> numByGeneration;
0500 numByGeneration.reserve(10);
0501
0502 for (auto it = begin; it != end; ++it) {
0503 auto& particle = *it;
0504 const auto pid = particle.particleId();
0505 if (pid.generation() >= numByGeneration.size()) {
0506 numByGeneration.resize(pid.generation() + 1, 0);
0507 }
0508 unsigned int nextSubParticle = numByGeneration[pid.generation()]++;
0509
0510 auto newPid = particle.particleId().setSubParticle(nextSubParticle);
0511 particle.setParticleId(newPid);
0512 }
0513 }
0514
0515 }
