File indexing completed on 2025-12-16 09:23:52
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "ActsExamples/Io/EDM4hep/EDM4hepSimInputConverter.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Utilities/Logger.hpp"
0015 #include "Acts/Utilities/ScopedTimer.hpp"
0016 #include "Acts/Utilities/Zip.hpp"
0017 #include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp"
0018 #include "ActsExamples/EventData/GeometryContainers.hpp"
0019 #include "ActsExamples/EventData/SimHit.hpp"
0020 #include "ActsExamples/EventData/SimParticle.hpp"
0021 #include "ActsExamples/Framework/ProcessCode.hpp"
0022 #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp"
0023 #include "ActsExamples/Utilities/Paths.hpp"
0024 #include "ActsFatras/EventData/Barcode.hpp"
0025 #include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp"
0026 #include "ActsPlugins/EDM4hep/EDM4hepUtil.hpp"
0027
0028 #include <algorithm>
0029 #include <cstdint>
0030 #include <iomanip>
0031 #include <map>
0032 #include <numeric>
0033 #include <ranges>
0034 #include <stdexcept>
0035
0036 #include <DD4hep/Detector.h>
0037 #include <boost/histogram.hpp>
0038 #include <boost/histogram/ostream.hpp>
0039 #include <edm4hep/MCParticle.h>
0040 #include <edm4hep/MCParticleCollection.h>
0041 #include <edm4hep/SimTrackerHit.h>
0042 #include <edm4hep/SimTrackerHitCollection.h>
0043 #include <podio/Frame.h>
0044 #include <podio/ObjectID.h>
0045
0046 namespace bh = boost::histogram;
0047
0048 namespace ActsExamples {
0049 namespace detail {
0050 struct ParticleInfo {
0051 std::size_t particleIndex;
0052
0053 };
0054
0055 }
0056
0057 EDM4hepSimInputConverter::EDM4hepSimInputConverter(const Config& config,
0058 Acts::Logging::Level level)
0059 : PodioInputConverter("EDM4hepSimInputConverter", level, config.inputFrame),
0060 m_cfg(config) {
0061 if (m_cfg.outputParticlesGenerator.empty()) {
0062 throw std::invalid_argument(
0063 "Missing output collection generator particles");
0064 }
0065 if (m_cfg.outputParticlesSimulation.empty()) {
0066 throw std::invalid_argument(
0067 "Missing output collection simulated particles");
0068 }
0069 if (m_cfg.outputSimHits.empty()) {
0070 throw std::invalid_argument("Missing output collection sim hits");
0071 }
0072
0073 if (m_cfg.outputSimVertices.empty()) {
0074 throw std::invalid_argument("Missing output collection sim vertices");
0075 }
0076
0077 m_outputParticlesGenerator.initialize(m_cfg.outputParticlesGenerator);
0078 m_outputParticlesSimulation.initialize(m_cfg.outputParticlesSimulation);
0079 m_outputSimHits.initialize(m_cfg.outputSimHits);
0080 m_outputSimHitAssociation.maybeInitialize(m_cfg.outputSimHitAssociation);
0081 m_outputSimVertices.initialize(m_cfg.outputSimVertices);
0082
0083 ACTS_INFO("Configured EDM4hepSimInputConverter:");
0084 auto printCut = [](std::optional<double> opt) {
0085 if (opt.has_value()) {
0086 return std::to_string(opt.value());
0087 } else {
0088 return std::string{"<none>"};
0089 }
0090 };
0091 ACTS_INFO("- particle r: [" << printCut(m_cfg.particleRMin) << ", "
0092 << printCut(m_cfg.particleRMax) << "] mm");
0093 ACTS_INFO("- particle z: [" << printCut(m_cfg.particleZMin) << ", "
0094 << printCut(m_cfg.particleZMax) << "] mm");
0095
0096 m_cfg.trackingGeometry->visitSurfaces([&](const auto* surface) {
0097 const auto* detElement =
0098 dynamic_cast<const ActsPlugins::DD4hepDetectorElement*>(
0099 surface->associatedDetectorElement());
0100
0101 if (detElement == nullptr) {
0102 ACTS_ERROR("Surface has no associated detector element");
0103 return;
0104 }
0105
0106 const auto translation = detElement->sourceElement()
0107 .nominal()
0108 .worldTransformation()
0109 .GetTranslation();
0110 Acts::Vector3 position;
0111 position << translation[0], translation[1], translation[2];
0112 position *= Acts::UnitConstants::cm;
0113
0114 m_surfaceMap.insert({detElement->sourceElement().key(), surface});
0115 });
0116 }
0117
0118 bool EDM4hepSimInputConverter::acceptParticle(
0119 const edm4hep::MCParticle& particle) const {
0120 double z = particle.getVertex().z * Acts::UnitConstants::mm;
0121 if (m_cfg.particleZMin.has_value() && z < m_cfg.particleZMin.value()) {
0122 ACTS_VERBOSE("Skipping particle with z=" << z << " < zMin="
0123 << m_cfg.particleZMin.value());
0124 return false;
0125 }
0126
0127 if (m_cfg.particleZMax.has_value() && z > m_cfg.particleZMax.value()) {
0128 ACTS_VERBOSE("Skipping particle with z=" << z << " > zMax="
0129 << m_cfg.particleZMax.value());
0130 return false;
0131 }
0132
0133 double r = std::hypot(particle.getVertex()[0], particle.getVertex()[1]) *
0134 Acts::UnitConstants::mm;
0135
0136 if (m_cfg.particleRMin.has_value() && r < m_cfg.particleRMin.value()) {
0137 ACTS_VERBOSE("Skipping particle with r=" << r << " < rMin="
0138 << m_cfg.particleRMin.value());
0139 return false;
0140 }
0141
0142 if (m_cfg.particleRMax.has_value() && r > m_cfg.particleRMax.value()) {
0143 ACTS_VERBOSE("Skipping particle with r=" << r << " > rMax="
0144 << m_cfg.particleRMax.value());
0145 return false;
0146 }
0147
0148 double pt = std::hypot(particle.getMomentum()[0], particle.getMomentum()[1]);
0149
0150 if (m_cfg.particlePtMin.has_value() && pt < m_cfg.particlePtMin.value()) {
0151 ACTS_VERBOSE("Skipping particle with pt=" << pt << " < ptMin="
0152 << m_cfg.particlePtMin.value());
0153 return false;
0154 }
0155
0156 if (m_cfg.particlePtMax.has_value() && pt > m_cfg.particlePtMax.value()) {
0157 ACTS_VERBOSE("Skipping particle with pt=" << pt << " > ptMax="
0158 << m_cfg.particlePtMax.value());
0159 return false;
0160 }
0161
0162 return true;
0163 }
0164
0165 bool EDM4hepSimInputConverter::particleOrDescendantsHaveHits(
0166 const edm4hep::MCParticle& particle,
0167 const std::function<std::uint8_t(const edm4hep::MCParticle&)>& getNumHits)
0168 const {
0169
0170 if (getNumHits(particle) > 0) {
0171 return true;
0172 }
0173
0174
0175 for (const auto& daughter : particle.getDaughters()) {
0176 if (particleOrDescendantsHaveHits(daughter, getNumHits)) {
0177 return true;
0178 }
0179 }
0180
0181 return false;
0182 }
0183
0184 namespace {
0185 bool isGeneratorStable(const edm4hep::MCParticle& particle) {
0186
0187 constexpr int kUndecayedPhysicalParticleStatus = 1;
0188 constexpr int kDecayedPhysicalParticleStatus = 2;
0189 int status = particle.getGeneratorStatus();
0190 return status == kUndecayedPhysicalParticleStatus ||
0191 status == kDecayedPhysicalParticleStatus;
0192 }
0193
0194 void findGeneratorStableParticles(
0195 const edm4hep::MCParticle& particle,
0196 std::vector<edm4hep::MCParticle>& outputParticles) {
0197
0198 if (particle.isCreatedInSimulation()) {
0199 return;
0200 }
0201
0202 if (isGeneratorStable(particle)) {
0203
0204 if (std::ranges::find(outputParticles, particle) == outputParticles.end()) {
0205 outputParticles.push_back(particle);
0206 }
0207 return;
0208 }
0209
0210 for (const auto& daughter : particle.getDaughters()) {
0211 findGeneratorStableParticles(daughter, outputParticles);
0212 }
0213 }
0214 }
0215
0216 ProcessCode EDM4hepSimInputConverter::convert(const AlgorithmContext& ctx,
0217 const podio::Frame& frame) const {
0218 ACTS_DEBUG("Reading EDM4hep inputs");
0219
0220 const auto& mcParticleCollection =
0221 frame.get<edm4hep::MCParticleCollection>(m_cfg.inputParticles);
0222
0223 ACTS_DEBUG("Total input particles: " << mcParticleCollection.size()
0224 << " particles");
0225
0226 std::vector<SimBarcode> unorderedParticlesInitial;
0227
0228
0229
0230
0231 std::vector<std::pair<Acts::Vector3, std::vector<int>>> primaryVertices;
0232
0233 std::size_t nVertexParticles = 0;
0234
0235 {
0236 Acts::ScopedTimer timer("Finding primary vertices", logger(),
0237 Acts::Logging::DEBUG);
0238 for (const auto& particle : mcParticleCollection) {
0239 if (particle.parents_size() > 0) {
0240
0241 continue;
0242 }
0243 const auto& vtx = particle.getEndpoint();
0244
0245
0246 Acts::Vector3 vtxPos = {vtx[0], vtx[1], vtx[2]};
0247 vtxPos *= Acts::UnitConstants::mm;
0248
0249
0250 auto it = std::ranges::find_if(primaryVertices, [&vtxPos](const auto& v) {
0251 return v.first == vtxPos;
0252 });
0253
0254 std::vector<int>* vertexParticles = nullptr;
0255
0256 if (it == primaryVertices.end()) {
0257 ACTS_VERBOSE("Found primary vertex at " << vtx.x << ", " << vtx.y
0258 << ", " << vtx.z);
0259 primaryVertices.emplace_back(vtxPos, std::vector<int>{});
0260 vertexParticles = &primaryVertices.back().second;
0261 } else {
0262 vertexParticles = &it->second;
0263 }
0264
0265 assert(vertexParticles != nullptr);
0266
0267 if (isGeneratorStable(particle)) {
0268 vertexParticles->push_back(particle.getObjectID().index);
0269 nVertexParticles += 1;
0270 }
0271
0272 nVertexParticles += particle.getDaughters().size();
0273 std::ranges::copy(
0274 particle.getDaughters() | std::views::transform([](const auto& p) {
0275 return p.getObjectID().index;
0276 }),
0277 std::back_inserter(*vertexParticles));
0278 }
0279 }
0280
0281 ACTS_DEBUG("Found " << primaryVertices.size() << " primary vertices with "
0282 << nVertexParticles << " outgoing particles in total");
0283
0284
0285 ParentRelationship parentRelationship;
0286
0287
0288
0289 std::unordered_map<int, detail::ParticleInfo> edm4hepParticleMap;
0290
0291 std::vector<std::uint16_t> numSimHits;
0292 numSimHits.resize(mcParticleCollection.size());
0293
0294 std::size_t nGeneratorParticles = 0;
0295 for (const auto& particle : mcParticleCollection) {
0296 if (!particle.isCreatedInSimulation()) {
0297 nGeneratorParticles += 1;
0298 }
0299 }
0300
0301 std::vector<const edm4hep::SimTrackerHitCollection*> simHitCollections;
0302 for (const auto& name : m_cfg.inputSimHits) {
0303 simHitCollections.push_back(
0304 &frame.get<edm4hep::SimTrackerHitCollection>(name));
0305 }
0306
0307
0308 for (const auto* inputHits : simHitCollections) {
0309 for (const auto& hit : *inputHits) {
0310 auto particle = ActsPlugins::EDM4hepUtil::getParticle(hit);
0311
0312 std::size_t index = particle.getObjectID().index;
0313
0314 auto& num = numSimHits.at(index);
0315 constexpr unsigned int maxNum =
0316 (1 << (sizeof(decltype(numSimHits)::value_type) * 8)) - 1;
0317
0318 if (num == maxNum) {
0319 throw std::runtime_error{"Hit count " + std::to_string(num) +
0320 " is at the limit of " +
0321 std::to_string(maxNum)};
0322 }
0323
0324 num += 1;
0325 }
0326 }
0327
0328 std::function<std::uint16_t(const edm4hep::MCParticle&)> getNumHits =
0329 [&numSimHits](const edm4hep::MCParticle& p) {
0330 return numSimHits.at(p.getObjectID().index);
0331 };
0332
0333 std::optional particlesGeneratorUnordered = std::vector<SimParticle>{};
0334
0335 std::size_t nPrimaryVertices = 0;
0336
0337 {
0338 Acts::ScopedTimer timer("Walking particle tree", logger(),
0339 Acts::Logging::DEBUG);
0340
0341 Acts::AveragingScopedTimer timerFindStable(
0342 "Finding generator-stable particles", logger(), Acts::Logging::DEBUG);
0343
0344 std::vector<edm4hep::MCParticle> generatorStableParticles;
0345
0346 for (const auto& [vtxPos, particles] : primaryVertices) {
0347 nPrimaryVertices += 1;
0348 ACTS_VERBOSE("Walking particle tree for primary vertex at "
0349 << vtxPos.x() << ", " << vtxPos.y() << ", " << vtxPos.z());
0350 std::size_t nParticles = 0;
0351 std::size_t nSecondaryVertices = 0;
0352 std::size_t maxGen = 0;
0353
0354 auto startSize = unorderedParticlesInitial.size();
0355
0356
0357
0358 generatorStableParticles.clear();
0359
0360 ACTS_VERBOSE("Finding generator stable particles in "
0361 << particles.size()
0362 << " particles recorded for this primary vertex");
0363
0364 {
0365 auto s = timerFindStable.sample();
0366 for (const auto& index : particles) {
0367 const auto& inParticle = mcParticleCollection.at(index);
0368 findGeneratorStableParticles(inParticle, generatorStableParticles);
0369 }
0370 }
0371
0372 ACTS_VERBOSE(
0373 "Have " << generatorStableParticles.size()
0374 << " generator stable particles for this primary vertex");
0375
0376 particlesGeneratorUnordered->reserve(particlesGeneratorUnordered->size() +
0377 generatorStableParticles.size());
0378
0379 for (const auto& genParticle : generatorStableParticles) {
0380 SimParticle particle =
0381 EDM4hepUtil::readParticle(genParticle)
0382 .withParticleId(SimBarcode()
0383 .withParticle(nParticles)
0384 .withVertexPrimary(nPrimaryVertices));
0385 particlesGeneratorUnordered->push_back(particle);
0386 ACTS_VERBOSE("+ add GEN particle " << particle);
0387 ACTS_VERBOSE(" - at " << particle.position().transpose());
0388
0389 const auto pid = particle.particleId();
0390 unorderedParticlesInitial.push_back(particle.particleId());
0391 edm4hepParticleMap[genParticle.getObjectID().index] =
0392 detail::ParticleInfo{.particleIndex =
0393 unorderedParticlesInitial.size() - 1};
0394 processChildren(genParticle, pid, unorderedParticlesInitial,
0395 parentRelationship, edm4hepParticleMap,
0396 nSecondaryVertices, maxGen, getNumHits);
0397 }
0398
0399 ACTS_VERBOSE("Primary vertex complete, produced "
0400 << (unorderedParticlesInitial.size() - startSize)
0401 << " particles and " << nSecondaryVertices
0402 << " secondary vertices in " << maxGen << " generations");
0403 std::span particleSpan{
0404 std::next(unorderedParticlesInitial.begin(), startSize),
0405 unorderedParticlesInitial.end()};
0406 setSubParticleIds(particleSpan);
0407 }
0408 }
0409
0410 ACTS_DEBUG("Found " << nGeneratorParticles << " generator particles and "
0411 << unorderedParticlesInitial.size() - nGeneratorParticles
0412 << " particles from simulation");
0413 ACTS_DEBUG("Found " << unorderedParticlesInitial.size()
0414 << " particles in total");
0415
0416 if (unorderedParticlesInitial.empty()) {
0417 ACTS_WARNING(
0418 "No particles found in the input file, this will likely cause issues "
0419 "if sim hits are converted");
0420 }
0421
0422 std::optional particlesSimulatedUnordered = std::vector<SimParticle>{};
0423
0424 auto convertPosition4 = [&](const edm4hep::MCParticle& inParticle) {
0425 Acts::Vector4 vtxPos4 = {inParticle.getVertex()[0],
0426 inParticle.getVertex()[1],
0427 inParticle.getVertex()[2], inParticle.getTime()};
0428 vtxPos4.head<3>() *= Acts::UnitConstants::mm;
0429 vtxPos4[3] *= Acts::UnitConstants::ns;
0430 return vtxPos4;
0431 };
0432
0433 std::size_t nSkipped = 0;
0434 {
0435 Acts::ScopedTimer timer("Converting particles", logger(),
0436 Acts::Logging::DEBUG);
0437
0438 for (const auto& inParticle : mcParticleCollection) {
0439 ACTS_VERBOSE("Converting particle:\n" << inParticle);
0440
0441 auto particleIt = edm4hepParticleMap.find(inParticle.getObjectID().index);
0442 if (particleIt == edm4hepParticleMap.end()) {
0443 nSkipped += 1;
0444 continue;
0445 }
0446
0447 const auto& [index] = particleIt->second;
0448
0449 if (!acceptParticle(inParticle) && getNumHits(inParticle) == 0) {
0450
0451 ACTS_VERBOSE(" - skipping particle (no hits AND not accepted)");
0452 continue;
0453 }
0454
0455 const auto& pid = unorderedParticlesInitial.at(index);
0456
0457
0458
0459 SimParticle particleSimulated = EDM4hepUtil::readParticle(inParticle);
0460 particleSimulated.setParticleId(pid);
0461 ACTS_VERBOSE("Have converted particle: " << particleSimulated);
0462
0463
0464
0465
0466 float time = inParticle.getTime() * Acts::UnitConstants::ns;
0467 ACTS_VERBOSE("Particle has " << inParticle.getDaughters().size()
0468 << " daughters");
0469 for (const auto& daughter : inParticle.getDaughters()) {
0470 if (!daughter.vertexIsNotEndpointOfParent()) {
0471 Acts::Vector4 pos4 = convertPosition4(daughter);
0472 time = static_cast<float>(pos4[Acts::eFreeTime]);
0473
0474 break;
0475 }
0476 }
0477
0478 particleSimulated.finalState().setPosition4(
0479 inParticle.getEndpoint()[0] * Acts::UnitConstants::mm,
0480 inParticle.getEndpoint()[1] * Acts::UnitConstants::mm,
0481 inParticle.getEndpoint()[2] * Acts::UnitConstants::mm, time);
0482
0483 Acts::Vector3 momentumFinal = {inParticle.getMomentumAtEndpoint()[0],
0484 inParticle.getMomentumAtEndpoint()[1],
0485 inParticle.getMomentumAtEndpoint()[2]};
0486 particleSimulated.finalState().setDirection(momentumFinal.normalized());
0487 particleSimulated.finalState().setAbsoluteMomentum(momentumFinal.norm());
0488
0489 ACTS_VERBOSE(
0490 "- Updated particle initial -> final, position: "
0491 << particleSimulated.initialState().fourPosition().transpose()
0492 << " -> "
0493 << particleSimulated.finalState().fourPosition().transpose());
0494 ACTS_VERBOSE(
0495 " momentum: "
0496 << particleSimulated.initialState().fourMomentum().transpose()
0497 << " -> "
0498 << particleSimulated.finalState().fourMomentum().transpose());
0499
0500 particlesSimulatedUnordered->push_back(particleSimulated);
0501 }
0502 }
0503
0504 ACTS_DEBUG("Converted " << particlesGeneratorUnordered->size()
0505 << " generator particles");
0506 ACTS_DEBUG("Converted " << particlesSimulatedUnordered->size()
0507 << " simulated particles");
0508
0509 ACTS_DEBUG("Skipped particles: " << nSkipped << " (no hits or not accepted)");
0510
0511 std::ranges::sort(*particlesGeneratorUnordered, detail::CompareParticleId{});
0512 std::ranges::sort(*particlesSimulatedUnordered, detail::CompareParticleId{});
0513
0514 SimParticleContainer particlesGenerator{particlesGeneratorUnordered->begin(),
0515 particlesGeneratorUnordered->end()};
0516
0517 particlesGeneratorUnordered.reset();
0518
0519 SimParticleContainer particlesSimulated{particlesSimulatedUnordered->begin(),
0520 particlesSimulatedUnordered->end()};
0521
0522 particlesSimulatedUnordered.reset();
0523
0524 std::optional simHitsUnordered = std::vector<SimHit>{};
0525 ACTS_DEBUG("Reading sim hits from " << m_cfg.inputSimHits.size()
0526 << " sim hit collections");
0527 {
0528 Acts::ScopedTimer timer("Reading sim hits", logger(), Acts::Logging::DEBUG);
0529
0530 const auto& vm = m_cfg.dd4hepDetector->dd4hepDetector().volumeManager();
0531
0532 auto geometryMapper = [&](std::uint64_t cellId) {
0533 ACTS_VERBOSE("CellID: " << cellId);
0534
0535 const auto detElement = vm.lookupDetElement(cellId);
0536
0537 ACTS_VERBOSE(" -> detElement: " << detElement.name());
0538 ACTS_VERBOSE(" -> id: " << detElement.id());
0539 ACTS_VERBOSE(" -> key: " << detElement.key());
0540
0541 Acts::Vector3 position;
0542 position
0543 << detElement.nominal().worldTransformation().GetTranslation()[0],
0544 detElement.nominal().worldTransformation().GetTranslation()[1],
0545 detElement.nominal().worldTransformation().GetTranslation()[2];
0546 position *= Acts::UnitConstants::cm;
0547
0548 ACTS_VERBOSE(" -> detElement position: " << position.transpose());
0549
0550 auto it = m_surfaceMap.find(detElement.key());
0551 if (it == m_surfaceMap.end()) {
0552 ACTS_ERROR("Unable to find surface for detElement "
0553 << detElement.name() << " with cellId " << cellId);
0554 throw std::runtime_error("Unable to find surface for detElement");
0555 }
0556 const auto* surface = it->second;
0557 if (surface == nullptr) {
0558 ACTS_ERROR("Unable to find surface for detElement "
0559 << detElement.name() << " with cellId " << cellId);
0560 throw std::runtime_error("Unable to find surface for detElement");
0561 }
0562 ACTS_VERBOSE(" -> surface: " << surface->geometryId());
0563 return surface->geometryId();
0564 };
0565
0566 auto particleMapper = [&](const auto& inParticle) {
0567 auto it = edm4hepParticleMap.find(inParticle.getObjectID().index);
0568 if (it == edm4hepParticleMap.end()) {
0569 ACTS_ERROR(
0570 "SimHit has source particle that we did not see "
0571 "before, particle=\n"
0572 << inParticle);
0573
0574 throw std::runtime_error(
0575 "SimHit has source particle that we did not see before, "
0576 "but expect to be here");
0577 }
0578
0579 auto pid = unorderedParticlesInitial.at(it->second.particleIndex);
0580 return pid;
0581 };
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591 std::size_t totalSimHitCount = 0;
0592 std::int32_t simHitIndex = 0;
0593
0594 if (m_outputSimHitAssociation.isInitialized()) {
0595 totalSimHitCount = std::accumulate(
0596 simHitCollections.begin(), simHitCollections.end(), 0u,
0597 [&](auto sum, const auto* coll) { return sum + coll->size(); });
0598
0599 if (totalSimHitCount >
0600 std::numeric_limits<decltype(simHitIndex)>::max()) {
0601 ACTS_ERROR(
0602 "Due to the way the conversion uses a 32bit integer to store the "
0603 "edm4hep sim hit index, the total number of sim hits across all "
0604 "input collections must be <= "
0605 << std::numeric_limits<decltype(simHitIndex)>::max() << ", but is "
0606 << totalSimHitCount);
0607 throw std::runtime_error("Total sim hit count is too large");
0608 }
0609 }
0610
0611 for (const auto& [inputHits, name] :
0612 Acts::zip(simHitCollections, m_cfg.inputSimHits)) {
0613 ACTS_VERBOSE("SimHit collection " << name << " has " << inputHits->size()
0614 << " hits");
0615
0616 simHitsUnordered->reserve(simHitsUnordered->size() + inputHits->size());
0617
0618 for (const auto& hit : *inputHits) {
0619 auto simHit = EDM4hepUtil::readSimHit(hit, particleMapper,
0620 geometryMapper, simHitIndex);
0621
0622 ACTS_VERBOSE("Converted sim hit for truth particle: "
0623 << simHit.particleId() << " at "
0624 << simHit.fourPosition().transpose() << " with time "
0625 << simHit.time());
0626
0627
0628 if (auto itSim = particlesSimulated.find(simHit.particleId());
0629 itSim != particlesSimulated.end()) {
0630 ACTS_VERBOSE("Found associated simulated particle");
0631 itSim->finalState().setNumberOfHits(
0632 itSim->finalState().numberOfHits() + 1);
0633 } else if (auto itGen = particlesGenerator.find(simHit.particleId());
0634 itGen != particlesGenerator.end()) {
0635 ACTS_VERBOSE("Found associated generator particle");
0636 itGen->finalState().setNumberOfHits(
0637 itGen->finalState().numberOfHits() + 1);
0638 } else {
0639 const auto& ptcl = ActsPlugins::EDM4hepUtil::getParticle(hit);
0640 ACTS_ERROR("SimHit (r="
0641 << Acts::VectorHelpers::perp(simHit.position())
0642 << ", z=" << simHit.position()[Acts::eFreePos2]
0643 << ") has source particle that we did not see before:\n"
0644 << ptcl);
0645 double particleR =
0646 std::hypot(ptcl.getVertex()[0], ptcl.getVertex()[1]) *
0647 Acts::UnitConstants::mm;
0648 double particleZ = ptcl.getVertex()[2] * Acts::UnitConstants::mm;
0649
0650 double particleREnd =
0651 std::hypot(ptcl.getEndpoint()[0], ptcl.getEndpoint()[1]) *
0652 Acts::UnitConstants::mm;
0653 double particleZEnd = ptcl.getEndpoint()[2] * Acts::UnitConstants::mm;
0654
0655 ACTS_ERROR("Particle loc: " << particleR << ", " << particleZ
0656 << " -> " << particleREnd << ", "
0657 << particleZEnd);
0658 continue;
0659 }
0660
0661 simHitsUnordered->push_back(std::move(simHit));
0662 simHitIndex += 1;
0663 }
0664 }
0665 }
0666
0667 ACTS_DEBUG("Read " << simHitsUnordered->size() << " sim hits in total");
0668
0669 std::ranges::sort(*simHitsUnordered, detail::CompareGeometryId{});
0670
0671 SimHitContainer simHits{simHitsUnordered->begin(), simHitsUnordered->end()};
0672 simHitsUnordered.reset();
0673
0674
0675
0676
0677
0678 if (m_outputSimHitAssociation.isInitialized()) {
0679 ActsPlugins::EDM4hepUtil::SimHitAssociation simHitAssociation;
0680 simHitAssociation.reserve(simHits.size());
0681
0682
0683 Acts::ScopedTimer timer("Building sim hit association", logger(),
0684 Acts::Logging::DEBUG);
0685 for (const auto&& [indexInColl, hit] : Acts::enumerate(simHits)) {
0686 std::size_t index = hit.index();
0687
0688 for (const auto&& [name, coll] :
0689 Acts::zip(m_cfg.inputSimHits, simHitCollections)) {
0690 if (index >= coll->size()) {
0691 index -= coll->size();
0692 continue;
0693 }
0694
0695 ACTS_VERBOSE("Hit assoc int -> ext #" << indexInColl << " -> "
0696 << coll->at(index).id() << " "
0697 << hit.position().transpose());
0698
0699 simHitAssociation.add(indexInColl, coll->at(index));
0700
0701 break;
0702 }
0703 }
0704
0705 if (simHitAssociation.size() != simHits.size()) {
0706 ACTS_ERROR("Sim hit association size " << simHitAssociation.size()
0707 << " does not match sim hit size "
0708 << simHits.size());
0709 throw std::runtime_error("Sim hit association size mismatch");
0710 }
0711
0712 m_outputSimHitAssociation(ctx, std::move(simHitAssociation));
0713 }
0714
0715 if (m_cfg.sortSimHitsInTime) {
0716 Acts::ScopedTimer timer("Sorting sim hits in time", logger(),
0717 Acts::Logging::DEBUG);
0718 std::multimap<ActsFatras::Barcode, std::size_t> hitsByParticle;
0719
0720 for (std::size_t i = 0; i < simHits.size(); ++i) {
0721 hitsByParticle.insert({simHits.nth(i)->particleId(), i});
0722 }
0723
0724 for (auto it = hitsByParticle.begin(), end = hitsByParticle.end();
0725 it != end; it = hitsByParticle.upper_bound(it->first)) {
0726 ACTS_VERBOSE("Particle " << it->first << " has "
0727 << hitsByParticle.count(it->first) << " hits");
0728
0729 std::vector<std::size_t> hitIndices;
0730 hitIndices.reserve(hitsByParticle.count(it->first));
0731 for (auto hitIndex : makeRange(hitsByParticle.equal_range(it->first))) {
0732 hitIndices.push_back(hitIndex.second);
0733 }
0734
0735 if (logger().doPrint(Acts::Logging::VERBOSE)) {
0736 ACTS_VERBOSE("Before sorting:");
0737 for (const auto& hitIdx : hitIndices) {
0738 ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index()
0739 << " t=" << simHits.nth(hitIdx)->time());
0740 }
0741 }
0742
0743 std::ranges::sort(hitIndices, {}, [&simHits](std::size_t h) {
0744 return simHits.nth(h)->time();
0745 });
0746
0747 for (std::size_t i = 0; i < hitIndices.size(); ++i) {
0748 auto& hit = *simHits.nth(hitIndices[i]);
0749
0750 hit = {hit.geometryId(), hit.particleId(),
0751 hit.fourPosition(), hit.momentum4Before(),
0752 hit.momentum4After(), static_cast<std::int32_t>(i)};
0753 }
0754
0755 if (logger().doPrint(Acts::Logging::VERBOSE)) {
0756 ACTS_VERBOSE("After sorting:");
0757 for (const auto& hitIdx : hitIndices) {
0758 ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index()
0759 << " t=" << simHits.nth(hitIdx)->time());
0760 }
0761 }
0762 }
0763 } else {
0764
0765
0766 for (auto& hit : simHits) {
0767
0768 hit = {hit.geometryId(), hit.particleId(), hit.fourPosition(),
0769 hit.momentum4Before(), hit.momentum4After(), -1};
0770 }
0771 }
0772
0773 std::vector<SimVertex> simVerticesUnordered;
0774
0775 auto maybeAddVertex = [&](const Acts::Vector4& vtxPos4,
0776
0777 SimVertexBarcode vtxId) -> SimVertex& {
0778 auto getMinDistance = [&]() {
0779 std::stringstream sstr;
0780 auto closestIt = std::ranges::min_element(
0781 simVerticesUnordered, {}, [&vtxPos4](const auto& v) {
0782 return (v.position4.template head<3>() - vtxPos4.template head<3>())
0783 .norm();
0784 });
0785
0786 if (closestIt != simVerticesUnordered.end()) {
0787 sstr << (closestIt->position4.head<3>() - vtxPos4.head<3>()).norm();
0788 } else {
0789 sstr << "[NONE]";
0790 }
0791
0792 return sstr.str();
0793 };
0794
0795 auto vertexIt =
0796 std::ranges::find_if(simVerticesUnordered, [&](const auto& v) {
0797 return (v.position4.template head<3>() - vtxPos4.template head<3>())
0798 .norm() < Acts::UnitConstants::mm * 1e-3 &&
0799 v.id == vtxId;
0800 });
0801
0802 SimVertex* vertex = nullptr;
0803
0804 if (vertexIt == simVerticesUnordered.end()) {
0805 ACTS_VERBOSE("Adding new vertex: position="
0806 << vtxPos4.template head<3>().transpose() << " id=" << vtxId
0807 << " (closest existing=" << getMinDistance() << ")");
0808 vertex = &simVerticesUnordered.emplace_back(vtxId, vtxPos4);
0809 } else {
0810 vertex = &*vertexIt;
0811 ACTS_VERBOSE("Reusing existing vertex: position="
0812 << vtxPos4.template head<3>().transpose()
0813 << " id=" << vertex->id);
0814 }
0815
0816 assert(vertex != nullptr);
0817
0818 return *vertex;
0819 };
0820
0821 {
0822 Acts::ScopedTimer timer("Finding source vertices for particles", logger(),
0823 Acts::Logging::DEBUG);
0824
0825 std::size_t nParticlesWithHits = 0;
0826 for (const auto& particle : particlesSimulated) {
0827
0828
0829 if (particle.finalState().numberOfHits() == 0) {
0830
0831 continue;
0832 }
0833 nParticlesWithHits += 1;
0834
0835 SimVertex& vertex = maybeAddVertex(
0836 particle.fourPosition(), SimVertexBarcode{particle.particleId()});
0837 vertex.outgoing.insert(particle.particleId());
0838 }
0839
0840 ACTS_DEBUG("Made " << simVerticesUnordered.size() << " vertices from "
0841 << nParticlesWithHits
0842 << " simulated particles with hits");
0843 }
0844
0845 std::ranges::sort(simVerticesUnordered, detail::CompareVertexId{});
0846
0847 ACTS_DEBUG("Converted number of vertices: " << simVerticesUnordered.size());
0848 SimVertexContainer simVertices{simVerticesUnordered.begin(),
0849 simVerticesUnordered.end()};
0850
0851 m_outputParticlesGenerator(ctx, std::move(particlesGenerator));
0852 m_outputParticlesSimulation(ctx, std::move(particlesSimulated));
0853 m_outputSimVertices(ctx, std::move(simVertices));
0854
0855 for (const auto&& [i, hit] : Acts::enumerate(simHits)) {
0856 ACTS_VERBOSE("- " << i << " " << hit.fourPosition().transpose());
0857 }
0858
0859 m_outputSimHits(ctx, std::move(simHits));
0860
0861 return ProcessCode::SUCCESS;
0862 }
0863
0864 void EDM4hepSimInputConverter::processChildren(
0865 const edm4hep::MCParticle& inParticle, SimBarcode parentId,
0866 std::vector<SimBarcode>& particles, ParentRelationship& parentRelationship,
0867 std::unordered_map<int, detail::ParticleInfo>& particleMap,
0868 std::size_t& nSecondaryVertices, std::size_t& maxGen,
0869 const std::function<std::uint8_t(const edm4hep::MCParticle&)>& getNumHits)
0870 const {
0871 constexpr auto indent = [&](std::size_t n) {
0872 std::string result;
0873 for (std::size_t i = 0; i < n; ++i) {
0874 result += " ";
0875 }
0876 return result;
0877 };
0878
0879 const std::size_t gen = parentId.generation();
0880 maxGen = std::max(maxGen, gen);
0881
0882 ACTS_VERBOSE(indent(gen) << " - processing daughters for input particle "
0883 << inParticle.id());
0884 ACTS_VERBOSE(indent(gen) << " -> found " << inParticle.daughters_size()
0885 << " daughter(s)");
0886
0887 bool parentDecayed =
0888 std::any_of(inParticle.daughters_begin(), inParticle.daughters_end(),
0889 [](const edm4hep::MCParticle& daughter) {
0890 return !daughter.vertexIsNotEndpointOfParent();
0891 });
0892 std::size_t secondaryVertex = 0;
0893
0894 Acts::Vector3 start{inParticle.getVertex().x, inParticle.getVertex().y,
0895 inParticle.getVertex().z};
0896 Acts::Vector3 end{inParticle.getEndpoint().x, inParticle.getEndpoint().y,
0897 inParticle.getEndpoint().z};
0898 double distance = (end - start).norm() * Acts::UnitConstants::mm;
0899
0900 constexpr double tolerance = Acts::s_onSurfaceTolerance;
0901 if (parentDecayed && distance > tolerance) {
0902 ACTS_VERBOSE(indent(gen) << " -> parent decays");
0903 secondaryVertex = ++nSecondaryVertices;
0904 }
0905
0906 std::size_t parentIndex = particles.size() - 1;
0907
0908 std::size_t nParticles = 0;
0909 for (const auto& daughter : inParticle.getDaughters()) {
0910 if (auto pIt = particleMap.find(daughter.getObjectID().index);
0911 pIt != particleMap.end()) {
0912 ACTS_WARNING("Skipping particle #"
0913 << daughter.getObjectID()
0914 << " that we've already processed, this should not happen "
0915 "at this stage.");
0916 continue;
0917 }
0918
0919
0920
0921 if (!acceptParticle(daughter) &&
0922 !particleOrDescendantsHaveHits(daughter, getNumHits)) {
0923 ACTS_VERBOSE(indent(gen + 1) << " - skipping particle " << daughter.id()
0924 << " (no hits AND not accepted)");
0925
0926 continue;
0927 }
0928
0929 SimParticle particle = EDM4hepUtil::readParticle(daughter);
0930
0931 auto pid = parentId.makeDescendant(nParticles);
0932 nParticles += 1;
0933 if (daughter.vertexIsNotEndpointOfParent()) {
0934
0935 } else {
0936
0937 pid = pid.withVertexSecondary(secondaryVertex);
0938 }
0939 particle.setParticleId(pid);
0940
0941 ACTS_VERBOSE(indent(particle.particleId().generation())
0942 << "+ add particle " << particle << " ("
0943 << particle.particleId() << ") from #"
0944 << daughter.getObjectID());
0945 ACTS_VERBOSE(indent(particle.particleId().generation())
0946 << " - generation: " << particle.particleId().generation());
0947 ACTS_VERBOSE(indent(particle.particleId().generation())
0948 << " - at " << particle.fourPosition().transpose());
0949 ACTS_VERBOSE(indent(particle.particleId().generation())
0950 << " - createdInSim: "
0951 << daughter.isCreatedInSimulation());
0952 ACTS_VERBOSE(indent(particle.particleId().generation())
0953 << " - vertexIsNotEndpointOfParent: "
0954 << daughter.vertexIsNotEndpointOfParent());
0955 ACTS_VERBOSE(indent(particle.particleId().generation())
0956 << " - isStopped: " << daughter.isStopped());
0957 ACTS_VERBOSE(indent(particle.particleId().generation())
0958 << " - endpoint: " << daughter.getEndpoint().x << ", "
0959 << daughter.getEndpoint().y << ", "
0960 << daughter.getEndpoint().z);
0961
0962 particles.push_back(particle.particleId());
0963 particleMap[daughter.getObjectID().index] =
0964 detail::ParticleInfo{.particleIndex = particles.size() - 1};
0965 parentRelationship[particles.size() - 1] = parentIndex;
0966
0967 processChildren(daughter, pid, particles, parentRelationship, particleMap,
0968 nSecondaryVertices, maxGen, getNumHits);
0969 }
0970 }
0971
0972 void EDM4hepSimInputConverter::setSubParticleIds(
0973 std::span<SimBarcode> particles) {
0974 std::vector<std::size_t> numByGeneration;
0975 numByGeneration.reserve(10);
0976
0977 for (auto& particle : particles) {
0978 if (particle.generation() >= numByGeneration.size()) {
0979 numByGeneration.resize(particle.generation() + 1, 0);
0980 }
0981 unsigned long nextSubParticle = numByGeneration[particle.generation()];
0982 numByGeneration[particle.generation()] += 1;
0983
0984 particle = particle.withSubParticle(nextSubParticle);
0985 }
0986 }
0987
0988 }
