diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index f32a5a13..7b49e5d4 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -43,6 +43,21 @@ inline float time_of_flight(const T& pos) { return std::sqrt((pos[0] * pos[0]) + (pos[1] * pos[1]) + (pos[2] * pos[2])) / TMath::C() * 1e6; } +// Index of the copied background particle a relation should point at, or -1 when +// the relation has to be left unset. An unset relation is any negative index, both +// podio::ObjectID::untracked (-1) and podio::ObjectID::invalid (-2), which is why +// oldIndex is tested against 0 rather than against -1. std::map::operator[] cannot +// be used for the lookup either: it default-constructs a 0 for a missing key, which +// would silently attach the hit to the first particle of the background event +// instead of leaving the relation unset. +inline int mapped_particle_index(const std::map& oldToNewMap, int oldIndex) { + if (oldIndex < 0) { + return -1; + } + const auto it = oldToNewMap.find(oldIndex); + return it == oldToNewMap.end() ? -1 : it->second; +} + std::pair OverlayTiming::define_time_windows(const std::string& collection_name) const { try { return {m_timeWindows.value().at(collection_name)[0], m_timeWindows.value().at(collection_name)[1]}; @@ -142,11 +157,15 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } // Fix relations to point to the new particles for (size_t i = 0; i < particles.size(); ++i) { - for (const auto& parent : particles[i].getParents()) { - oparticles[i].addToParents(oparticles[parent.getObjectID().index]); + for (const auto& parent : particles.at(i).getParents()) { + if (const auto index = parent.getObjectID().index; index >= 0) { + oparticles.at(i).addToParents(oparticles.at(index)); + } } - for (const auto& daughter : particles[i].getDaughters()) { - oparticles[i].addToDaughters(oparticles[daughter.getObjectID().index]); + for (const auto& daughter : particles.at(i).getDaughters()) { + if (const auto index = daughter.getObjectID().index; index >= 0) { + oparticles.at(i).addToDaughters(oparticles.at(index)); + } } } @@ -160,8 +179,8 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, const float tof = time_of_flight(simTrackerHit.getPosition()); if ((simTrackerHit.getTime() > this_start + tof) && (simTrackerHit.getTime() < this_stop + tof)) { auto nhit = simTrackerHit.clone(false); - if (simTrackerHit.getParticle().getObjectID().index != -1) { - nhit.setParticle(oparticles[simTrackerHit.getParticle().getObjectID().index]); + if (const auto index = simTrackerHit.getParticle().getObjectID().index; index >= 0) { + nhit.setParticle(oparticles.at(index)); } ocoll.push_back(nhit); } @@ -175,6 +194,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, const auto& coll = simCaloHits[i]; const auto name = inputLocations(SIMCALOHIT_INDEX_POSITION)[i]; const auto [this_start, this_stop] = define_time_windows(name); + // operator[] on purpose: this is where the entry for this collection is created auto& calHitMap = cellIDsMap[i]; auto& caloHitContribs = ocaloHitContribs[i]; for (const auto&& simCaloHit : *coll) { @@ -188,14 +208,18 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, within_time_window = true; // TODO: Make sure a contribution is not added twice auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[contrib.getParticle().getObjectID().index]); + // The contribution may have no particle attached, in which case the + // relation is left unset rather than indexed with a negative index. + if (const auto index = contrib.getParticle().getObjectID().index; index >= 0) { + newContrib.setParticle(oparticles.at(index)); + } thisContribs.push_back(caloHitContribs.size()); caloHitContribs.push_back(std::move(newContrib)); } if (within_time_window) { auto newhit = simCaloHit.clone(false); for (const auto& contrib : thisContribs) { - newhit.addToContributions(caloHitContribs[contrib]); + newhit.addToContributions(caloHitContribs.at(contrib)); } calHitMap.emplace(simCaloHit.getCellID(), std::move(newhit)); } @@ -234,14 +258,14 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // Overlay the background events to each bunchcrossing in the bunch train for (int bxInTrain = 0; bxInTrain < m_NBunchTrain; ++bxInTrain) { - const int BX_number_in_train = permutation.at(bxInTrain); + const int BX_number_in_train = permutation[bxInTrain]; int NOverlay_to_this_BX = 0; - if (m_Poisson[groupIndex]) { - NOverlay_to_this_BX = std::poisson_distribution<>(m_Noverlay[groupIndex])(rng_engine); + if (m_Poisson.value()[groupIndex]) { + NOverlay_to_this_BX = std::poisson_distribution<>(m_Noverlay.value()[groupIndex])(rng_engine); } else { - NOverlay_to_this_BX = m_Noverlay[groupIndex]; + NOverlay_to_this_BX = m_Noverlay.value()[groupIndex]; } debug() << "Will overlay " << NOverlay_to_this_BX << " events to BX number " << BX_number_in_train + physBX @@ -275,15 +299,15 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, const auto& bgParticles = backgroundEvent.get(m_MCParticleCollectionName); int j = oparticles.size(); for (size_t i = 0; i < bgParticles.size(); ++i) { - auto npart = bgParticles[i].clone(false); + auto npart = bgParticles.at(i).clone(false); - npart.setTime(bgParticles[i].getTime() + timeOffset); + npart.setTime(bgParticles.at(i).getTime() + timeOffset); npart.setOverlay(true); oparticles.push_back(npart); - for (const auto& parent : bgParticles[i].getParents()) { + for (const auto& parent : bgParticles.at(i).getParents()) { parentDaughterMap[j].first.push_back(parent.getObjectID().index); } - for (const auto& daughter : bgParticles[i].getDaughters()) { + for (const auto& daughter : bgParticles.at(i).getDaughters()) { parentDaughterMap[j].second.push_back(daughter.getObjectID().index); } oldToNewMap[i] = j; @@ -292,19 +316,21 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, for (const auto& [index, parentsDaughters] : parentDaughterMap) { const auto& [parents, daughters] = parentsDaughters; for (const auto& parent : parents) { - if (parentDaughterMap.find(oldToNewMap[parent]) == parentDaughterMap.end()) { + const auto newIndex = mapped_particle_index(oldToNewMap, parent); + if (newIndex < 0 || parentDaughterMap.find(newIndex) == parentDaughterMap.end()) { // warning() << "Parent " << parent << " not found in background event" << endmsg; continue; } - oparticles[index].addToParents(oparticles[oldToNewMap[parent]]); + oparticles.at(index).addToParents(oparticles.at(newIndex)); } for (const auto& daughter : daughters) { - if (parentDaughterMap.find(oldToNewMap[daughter]) == parentDaughterMap.end()) { + const auto newIndex = mapped_particle_index(oldToNewMap, daughter); + if (newIndex < 0 || parentDaughterMap.find(newIndex) == parentDaughterMap.end()) { // warning() << "Parent " << daughter << " not found in background event" << endmsg; continue; } // info() << "Adding (daughter) " << daughter << " to " << index << endmsg; - oparticles[index].addToDaughters(oparticles[oldToNewMap[daughter]]); + oparticles.at(index).addToDaughters(oparticles.at(newIndex)); } } @@ -332,7 +358,10 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, auto nhit = simTrackerHit.clone(false); nhit.setOverlay(true); nhit.setTime(simTrackerHit.getTime() + timeOffset); - nhit.setParticle(oparticles[oldToNewMap[simTrackerHit.getParticle().getObjectID().index]]); + if (const auto index = mapped_particle_index(oldToNewMap, simTrackerHit.getParticle().getObjectID().index); + index >= 0) { + nhit.setParticle(oparticles.at(index)); + } ocoll.push_back(nhit); } } @@ -363,7 +392,10 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, add = true; // TODO: Make sure a contribution is not added twice auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + if (const auto index = mapped_particle_index(oldToNewMap, contrib.getParticle().getObjectID().index); + index >= 0) { + newContrib.setParticle(oparticles.at(index)); + } newContrib.setTime(contrib.getTime() + timeOffset); calhit.addToContributions(newContrib); calHitContribs.push_back(newContrib); @@ -382,7 +414,10 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { // TODO: Make sure a contribution is not added twice auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + if (const auto index = mapped_particle_index(oldToNewMap, contrib.getParticle().getObjectID().index); + index >= 0) { + newContrib.setParticle(oparticles.at(index)); + } newContrib.setTime(contrib.getTime() + timeOffset); calhit.addToContributions(newContrib); calHitContribs.push_back(newContrib); diff --git a/test/k4FWCoreTest/CMakeLists.txt b/test/k4FWCoreTest/CMakeLists.txt index 6f2e1dcc..0c6ddbd7 100644 --- a/test/k4FWCoreTest/CMakeLists.txt +++ b/test/k4FWCoreTest/CMakeLists.txt @@ -36,6 +36,7 @@ set(k4fwcoretest_plugin_sources src/components/ExampleFunctionalProducer.cpp src/components/ExampleFunctionalProducerMultiple.cpp src/components/ExampleFunctionalProducerRuntimeCollections.cpp + src/components/ExampleSimHitsWithoutParticles.cpp src/components/ExampleFunctionalTransformer.cpp src/components/ExampleFunctionalTransformerHist.cpp src/components/ExampleFunctionalTransformerMultiple.cpp @@ -249,7 +250,7 @@ add_test_fwcore(CheckLoadedFilesHaveCorrectDunderFile options/checkLoadedFilePro # appear in the output add_test_fwcore(CheckLoadedFileCorrectPathOnError options/checkLoadedFileProperties.py --with-error) set_tests_properties(CheckLoadedFileCorrectPathOnError - PROPERTIES PASS_REGULAR_EXPRESSION [=[ File ".*/test/k4FWCoreTest/options/checkLoadedFileProperties.py", line 36, in ]=] + PROPERTIES PASS_REGULAR_EXPRESSION [=[ File ".*/test/k4FWCoreTest/options/checkLoadedFileProperties.py", line 35, in ]=] ) add_test_fwcore(ParticleIDMetadataFramework options/ExampleParticleIDMetadata.py PROPERTIES FIXTURES_SETUP ParticleIDMetadataFile) @@ -274,6 +275,9 @@ set_tests_properties(InvalidOutputCommandsNoCrash PROPERTIES FIXTURES_REQUIRED P add_test_fwcore(OverlayTiming options/TestOverlayTiming.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_REQUIRED ProducerMultipleFile) +add_test_fwcore(SimHitsWithoutParticles options/ExampleSimHitsWithoutParticles.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_SETUP SimHitsWithoutParticlesFile) +add_test_fwcore(OverlayTimingNoParticles options/TestOverlayTimingNoParticles.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_REQUIRED SimHitsWithoutParticlesFile) + add_test(NAME check_broken_pipe COMMAND bash -c "[ $(${K4RUN} options/ExampleFunctionalProducer.py | head -n 2 | wc -l) = 2 ]" WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR} diff --git a/test/k4FWCoreTest/options/ExampleFunctionalProducerMultiple.py b/test/k4FWCoreTest/options/ExampleFunctionalProducerMultiple.py index 6d9099f0..1a4ea370 100644 --- a/test/k4FWCoreTest/options/ExampleFunctionalProducerMultiple.py +++ b/test/k4FWCoreTest/options/ExampleFunctionalProducerMultiple.py @@ -25,7 +25,6 @@ from k4FWCore import ApplicationMgr, IOSvc from Configurables import EventDataSvc - iosvc = IOSvc("IOSvc") iosvc.Output = "functional_producer_multiple.root" # Collections can be dropped diff --git a/test/k4FWCoreTest/options/ExampleSimHitsWithoutParticles.py b/test/k4FWCoreTest/options/ExampleSimHitsWithoutParticles.py new file mode 100644 index 00000000..33387d53 --- /dev/null +++ b/test/k4FWCoreTest/options/ExampleSimHitsWithoutParticles.py @@ -0,0 +1,38 @@ +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +# Writes sim hits whose MCParticle relations are unset, used as both signal and +# background by TestOverlayTimingNoParticles.py. + +from Gaudi.Configuration import INFO +from Configurables import EventDataSvc, ExampleSimHitsWithoutParticles +from k4FWCore import ApplicationMgr, IOSvc + +producer = ExampleSimHitsWithoutParticles("ExampleSimHitsWithoutParticles") + +iosvc = IOSvc("IOSvc") +iosvc.Output = "sim_hits_without_particles.root" + +ApplicationMgr( + TopAlg=[producer], + EvtSel="NONE", + EvtMax=3, + ExtSvc=[EventDataSvc("EventDataSvc")], + OutputLevel=INFO, +) diff --git a/test/k4FWCoreTest/options/TestAlgorithmWithTFile.py b/test/k4FWCoreTest/options/TestAlgorithmWithTFile.py index 55cdcb9b..3ebd3391 100644 --- a/test/k4FWCoreTest/options/TestAlgorithmWithTFile.py +++ b/test/k4FWCoreTest/options/TestAlgorithmWithTFile.py @@ -21,7 +21,6 @@ from k4FWCore import ApplicationMgr, IOSvc from Configurables import k4FWCoreTest_AlgorithmWithTFile, EventDataSvc - producer = k4FWCoreTest_AlgorithmWithTFile() iosvc = IOSvc() diff --git a/test/k4FWCoreTest/options/TestOverlayTimingNoParticles.py b/test/k4FWCoreTest/options/TestOverlayTimingNoParticles.py new file mode 100644 index 00000000..cc25fc47 --- /dev/null +++ b/test/k4FWCoreTest/options/TestOverlayTimingNoParticles.py @@ -0,0 +1,65 @@ +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +# Regression test: OverlayTiming must cope with SimTrackerHits and +# CaloHitContributions whose MCParticle relation is unset, both in the signal +# event and in the background event. Indexing the output particle collection +# with the -1 of an unset ObjectID used to produce a broken relation that +# crashed when the contributions were written out. + +from Gaudi.Configuration import INFO +from Configurables import EventDataSvc, EventHeaderCreator, OverlayTiming, UniqueIDGenSvc +from k4FWCore import ApplicationMgr, IOSvc + +uid_svc = UniqueIDGenSvc("UniqueIDGenSvc") + +iosvc = IOSvc("IOSvc") +iosvc.Input = "sim_hits_without_particles.root" +iosvc.Output = "overlay_no_particles_output.root" + +header = EventHeaderCreator("EventHeaderCreator") + +overlay = OverlayTiming("OverlayTiming") +overlay.MCParticles = "MCParticles" +overlay.SimTrackerHits = ["SimTrackerHits"] +overlay.SimCalorimeterHits = ["SimCalorimeterHits"] +overlay.OutputMCParticles = "OverlayMCParticles" +overlay.OutputSimTrackerHits = ["OverlaySimTrackerHits"] +overlay.OutputSimCalorimeterHits = ["OverlaySimCalorimeterHits"] +overlay.OutputCaloHitContributions = ["OverlayCaloHitContributions"] +overlay.BackgroundMCParticleCollectionName = "MCParticles" +# The same file is used as signal and as background, so the unset relations are +# exercised on both the signal-copy and the background-merge path. +overlay.BackgroundFileNames = [["sim_hits_without_particles.root"]] +overlay.NumberBackground = [1] +overlay.Poisson_random_NOverlay = [False] +overlay.NBunchtrain = 1 +overlay.AllowReusingBackgroundFiles = True +overlay.TimeWindows = { + "SimTrackerHits": [-10000, 10000], + "SimCalorimeterHits": [-10000, 10000], +} + +ApplicationMgr( + TopAlg=[header, overlay], + EvtSel="NONE", + EvtMax=3, + ExtSvc=[EventDataSvc("EventDataSvc"), uid_svc], + OutputLevel=INFO, +) diff --git a/test/k4FWCoreTest/options/checkExampleEventData.py b/test/k4FWCoreTest/options/checkExampleEventData.py index fccc6557..6ea69cd6 100644 --- a/test/k4FWCoreTest/options/checkExampleEventData.py +++ b/test/k4FWCoreTest/options/checkExampleEventData.py @@ -24,7 +24,6 @@ from Configurables import k4FWCoreTest_CheckExampleEventData from k4FWCore import ApplicationMgr, IOSvc - parser.add_argument( "--collections", action="extend", diff --git a/test/k4FWCoreTest/options/checkLoadedFileProperties.py b/test/k4FWCoreTest/options/checkLoadedFileProperties.py index 80a7589b..3c28aff1 100644 --- a/test/k4FWCoreTest/options/checkLoadedFileProperties.py +++ b/test/k4FWCoreTest/options/checkLoadedFileProperties.py @@ -22,7 +22,6 @@ from k4FWCore.parseArgs import parser - parser.add_argument( "--with-error", action="store_true", diff --git a/test/k4FWCoreTest/src/components/ExampleSimHitsWithoutParticles.cpp b/test/k4FWCoreTest/src/components/ExampleSimHitsWithoutParticles.cpp new file mode 100644 index 00000000..65e91675 --- /dev/null +++ b/test/k4FWCoreTest/src/components/ExampleSimHitsWithoutParticles.cpp @@ -0,0 +1,83 @@ +/* + * Copyright (c) 2014-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "k4FWCore/Producer.h" + +#include "edm4hep/CaloHitContributionCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" + +#include +#include + +// Produces sim hits whose MCParticle relations are deliberately left unset. This +// is what the hits of an event that has already been through an overlay look like +// when the background particles were not kept, and it is also simply what a +// SimCalorimeterHit built without a particle looks like. OverlayTiming has to +// leave such relations unset instead of indexing its particle collection with the +// negative index of the unset ObjectID. + +using retType = std::tuple; + +struct ExampleSimHitsWithoutParticles final : k4FWCore::Producer { + ExampleSimHitsWithoutParticles(const std::string& name, ISvcLocator* svcLoc) + : Producer(name, svcLoc, {}, + {KeyValues("OutputCollectionParticles", {"MCParticles"}), + KeyValues("OutputCollectionSimTrackerHits", {"SimTrackerHits"}), + KeyValues("OutputCollectionSimCalorimeterHits", {"SimCalorimeterHits"}), + KeyValues("OutputCollectionCaloHitContributions", {"CaloHitContributions"})}) {} + + retType operator()() const override { + // A couple of particles, so that the output particle collection is not empty + // and a wrong index would go unnoticed. + auto particles = edm4hep::MCParticleCollection(); + const edm4hep::Vector3d v{0, 0, 0}; + particles.create(1, 2, 3, 4.f, 5.f, 6.f, v, v, v); + particles.create(2, 3, 4, 5.f, 6.f, 7.f); + + auto simTrackerHits = edm4hep::SimTrackerHitCollection(); + auto simCaloHits = edm4hep::SimCalorimeterHitCollection(); + auto contributions = edm4hep::CaloHitContributionCollection(); + + for (int i = 0; i < 3; ++i) { + // No setParticle call: the relation stays unset. + auto trackerHit = simTrackerHits.create(); + trackerHit.setCellID(i + 1); + trackerHit.setTime(1.f); + trackerHit.setPosition({0, 0, 0}); + + auto contribution = contributions.create(); + contribution.setEnergy(1.f); + contribution.setTime(1.f); + + auto caloHit = simCaloHits.create(); + caloHit.setCellID(i + 1); + caloHit.setEnergy(1.f); + caloHit.setPosition({0, 0, 0}); + caloHit.addToContributions(contribution); + } + + return std::make_tuple(std::move(particles), std::move(simTrackerHits), std::move(simCaloHits), + std::move(contributions)); + } +}; + +DECLARE_COMPONENT(ExampleSimHitsWithoutParticles)