diff --git a/k4EDM4hep2LcioConv/CMakeLists.txt b/k4EDM4hep2LcioConv/CMakeLists.txt index 7495cf52..6c923574 100644 --- a/k4EDM4hep2LcioConv/CMakeLists.txt +++ b/k4EDM4hep2LcioConv/CMakeLists.txt @@ -10,7 +10,9 @@ target_include_directories(k4EDM4hep2LcioConv PUBLIC target_link_libraries(k4EDM4hep2LcioConv PUBLIC LCIO::lcio - EDM4HEP::edm4hep) + EDM4HEP::edm4hep + EDM4HEP::utils +) set(public_headers include/${PROJECT_NAME}/k4EDM4hep2LcioConv.h diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h index eb3b3b75..ae271a4a 100644 --- a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h +++ b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h @@ -33,6 +33,7 @@ namespace edm4hep { #endif #include #include +#include #include "podio/Frame.h" @@ -59,6 +60,7 @@ namespace edm4hep { #include #include +#include // Preprocessor symbol that can be used in downstream code to switch on the // namespace for the conversion @@ -85,8 +87,45 @@ namespace EDM4hep2LCIOConv { ObjectMapT vertices {}; ObjectMapT recoParticles {}; ObjectMapT mcParticles {}; + ObjectMapT particleIDs {}; }; + /// The minimal necessary information to do late conversions of ParticleIDs, + /// which is necessary to have consistent algorithmType information in the + /// conversion + struct ParticleIDConvData { + std::string name; + const edm4hep::ParticleIDCollection* coll; + std::optional metadata; + }; + + /// Sort the ParticleIDs according to their algorithmType. + /// + /// This sorting allows for a fully consistent roundtrip conversion under the + /// following conditions: + /// - All ParticleIDs are converted (where all means all that belong to a + /// given ReconstructedParticle collection) + /// - All of these conversions had the necessary metadata available + /// + /// In case these conditions are not met, the assigned algorithmTypes might + /// differ between LCIO and EDM4hep, but the metadata that is set will be + /// consistent for usage. + void sortParticleIDs(std::vector& pidCollections); + + /// Attach the meta information for one ParticleID collection to the passed + /// LCIO Event + /// + /// This returns the algorithmID according to the PIDHandler that is used to + /// attach the information to the LCIO reconstructed particle collection. If + /// there is no fitting reconstructed particle collection to attach this to + /// but the meta data in the ParticleIDConvData is valid, the algoType of that + /// will be returned. If there is no reconstructed particle collection or the + /// meta data is invalid an empty optional is returned. + std::optional attachParticleIDMetaData( + IMPL::LCEventImpl* lcEvent, + const podio::Frame& edmEvent, + const ParticleIDConvData& pidCollMetaInfo); + /** * Convert EDM4hep Tracks to LCIO. Simultaneously populate the mapping from * EDM4hep to LCIO objects for relation resolving in a second step. @@ -321,6 +360,15 @@ namespace EDM4hep2LCIOConv { return convertMCParticles(mcparticle_coll, mc_particles_vec).release(); } + /** + * Convert EDM4hep ParticleIDs to LCIO. NOTE: Since ParticleIDs cannot live in + * their own collections in LCIO this simply populates the pidMap that maps + * LCIO to EDM4hep particleIDs. **This just converts the data it is crucial to + * also resolve the relations afterwards!** + */ + template + void convertParticleIDs(const edm4hep::ParticleIDCollection* const edmCollection, PidMapT& pidMap, const int algoId); + /** * Convert EDM4hep EventHeader to LCIO. This will directly populate the * corresponding information in the passed LCEvent. The input collection needs @@ -389,6 +437,9 @@ namespace EDM4hep2LCIOConv { template void resolveRelationsClusters(ClusterMapT& clustersMap, const CaloHitMapT& caloHitMap); + template + void resolveRelationsParticleIDs(PidMapT& pidMap, const RecoParticleMapT& recoMap); + /** * Resolve all relations in all converted objects that are held in the map. * Dispatch to the correpsonding implementation for all the types that have diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp index fbdadc5b..fabae290 100644 --- a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp +++ b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp @@ -10,7 +10,6 @@ namespace EDM4hep2LCIOConv { TrackMapT& trackMap) { auto tracks = std::make_unique(lcio::LCIO::TRACK); - // Loop over EDM4hep tracks converting them to lcio tracks. for (const auto& edm_tr : (*edmCollection)) { if (edm_tr.isAvailable()) { @@ -392,23 +391,6 @@ namespace EDM4hep2LCIOConv { subdetEnergies.push_back(edmEnergy); } - // Convert ParticleIDs associated to the recoparticle - for (const auto& edm_pid : edm_cluster.getParticleIDs()) { - if (edm_pid.isAvailable()) { - auto* lcio_pid = new lcio::ParticleIDImpl; - - lcio_pid->setType(edm_pid.getType()); - lcio_pid->setPDG(edm_pid.getPDG()); - lcio_pid->setLikelihood(edm_pid.getLikelihood()); - lcio_pid->setAlgorithmType(edm_pid.getAlgorithmType()); - for (const auto& param : edm_pid.getParameters()) { - lcio_pid->addParameter(param); - } - - lcio_cluster->addParticleID(lcio_pid); - } - } - // Add LCIO and EDM4hep pair collections to vec k4EDM4hep2LcioConv::detail::mapInsert(lcio_cluster, edm_cluster, clusterMap); @@ -473,43 +455,6 @@ namespace EDM4hep2LCIOConv { float rp[3] = {edm_rp.getReferencePoint()[0], edm_rp.getReferencePoint()[1], edm_rp.getReferencePoint()[2]}; lcio_recp->setReferencePoint(rp); lcio_recp->setGoodnessOfPID(edm_rp.getGoodnessOfPID()); - - // Convert ParticleIDs associated to the recoparticle - for (const auto& edm_pid : edm_rp.getParticleIDs()) { - if (edm_pid.isAvailable()) { - auto* lcio_pid = new lcio::ParticleIDImpl; - - lcio_pid->setType(edm_pid.getType()); - lcio_pid->setPDG(edm_pid.getPDG()); - lcio_pid->setLikelihood(edm_pid.getLikelihood()); - lcio_pid->setAlgorithmType(edm_pid.getAlgorithmType()); - for (const auto& param : edm_pid.getParameters()) { - lcio_pid->addParameter(param); - } - - lcio_recp->addParticleID(lcio_pid); - } - } - - // Link sinlge associated Particle - auto edm_pid_used = edm_rp.getParticleIDUsed(); - if (edm_pid_used.isAvailable()) { - for (const auto& lcio_pid : lcio_recp->getParticleIDs()) { - bool is_same = true; - is_same = is_same && (lcio_pid->getType() == edm_pid_used.getType()); - is_same = is_same && (lcio_pid->getPDG() == edm_pid_used.getPDG()); - is_same = is_same && (lcio_pid->getLikelihood() == edm_pid_used.getLikelihood()); - is_same = is_same && (lcio_pid->getAlgorithmType() == edm_pid_used.getAlgorithmType()); - for (auto i = 0u; i < edm_pid_used.parameters_size(); ++i) { - is_same = is_same && (edm_pid_used.getParameters(i) == lcio_pid->getParameters()[i]); - } - if (is_same) { - lcio_recp->setParticleIDUsed(lcio_pid); - break; - } - } - } - // Add LCIO and EDM4hep pair collections to vec k4EDM4hep2LcioConv::detail::mapInsert(lcio_recp, edm_rp, recoparticles_vec); @@ -574,6 +519,22 @@ namespace EDM4hep2LCIOConv { return mcparticles; } + template + void convertParticleIDs(const edm4hep::ParticleIDCollection* const edmCollection, PidMapT& pidMap, const int algoId) + { + for (const auto& edmPid : (*edmCollection)) { + auto [lcioPid, _] = k4EDM4hep2LcioConv::detail::mapInsert(new lcio::ParticleIDImpl(), edmPid, pidMap).first; + + lcioPid->setType(edmPid.getType()); + lcioPid->setPDG(edmPid.getPDG()); + lcioPid->setLikelihood(edmPid.getLikelihood()); + lcioPid->setAlgorithmType(algoId); + for (const auto& param : edmPid.getParameters()) { + lcioPid->addParameter(param); + } + } + } + template void resolveRelationsMCParticles(MCParticleMapT& mcParticlesMap, const MCParticleLookupMapT& lookupMap) { @@ -762,6 +723,21 @@ namespace EDM4hep2LCIOConv { } } + template + void resolveRelationsParticleIDs(PidMapT& pidMap, const RecoParticleMapT& recoMap) + { + for (auto& [lcioPid, edmPid] : pidMap) { + const auto edmReco = edmPid.getParticle(); + const auto lcioReco = k4EDM4hep2LcioConv::detail::mapLookupFrom(edmReco, recoMap); + if (lcioReco) { + lcioReco.value()->addParticleID(lcioPid); + } + else { + std::cerr << "Cannot find a reconstructed particle to attach a ParticleID to" << std::endl; + } + } + } + template void resolveRelations(ObjectMappingT& collection_pairs) { @@ -783,6 +759,7 @@ namespace EDM4hep2LCIOConv { lookup_pairs.vertices, lookup_pairs.clusters, lookup_pairs.tracks); + resolveRelationsParticleIDs(lookup_pairs.particleIDs, update_pairs.recoParticles); resolveRelationsVertices(update_pairs.vertices, lookup_pairs.recoParticles); resolveRelationsSimCaloHit(update_pairs.simCaloHits, lookup_pairs.mcParticles); resolveRelationsSimTrackerHits(update_pairs.simTrackerHits, lookup_pairs.mcParticles); diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.h b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.h index 8a415536..2b0272ef 100644 --- a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.h +++ b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.h @@ -35,6 +35,7 @@ namespace edm4hep { #endif #include "edm4hep/TrackerHitPlaneCollection.h" #include "edm4hep/VertexCollection.h" +#include "edm4hep/utils/ParticleIDUtils.h" // LCIO #include @@ -90,7 +91,6 @@ namespace LCIO2EDM4hepConv { ObjectMapT recoParticles {}; ObjectMapT mcParticles {}; ObjectMapT trackerHitPlanes {}; - ObjectMapT particleIDs {}; }; using CollNamePair = std::tuple>; @@ -172,6 +172,21 @@ namespace LCIO2EDM4hepConv { inline edm4hep::Vector3f Vector3fFrom(const EVENT::FloatVec& v) { return edm4hep::Vector3f(v[0], v[1], v[2]); } + /** + * Get the name of a ParticleID collection from the name of the reco + * collection (from which it is created) and the PID algorithm name. + */ + inline std::string getPIDCollName(const std::string& recoCollName, const std::string& algoName) + { + return recoCollName + "_PID_" + algoName; + } + + /** + * Get the meta information for all particle id collections that are available + * from the PIDHandler + */ + std::vector getPIDMetaInfo(const EVENT::LCCollection* recoColl); + /** * Convert a TrackState */ @@ -200,16 +215,13 @@ namespace LCIO2EDM4hepConv { * Convert a ReconstructedParticle collection and return the resulting collection. * Simultaneously populates the mapping from LCIO to EDM4hep objects. * - * NOTE: Also populates a ParticleID collection, as those are persisted as + * NOTE: Also populates ParticleID collections, as those are persisted as * part of the ReconstructedParticles in LCIO. The name of this collection is - * _particleIDs + * _PID_ (see getPIDCollName) */ - template - std::vector convertReconstructedParticles( - const std::string& name, - EVENT::LCCollection* LCCollection, - RecoMapT& recoparticlesMap, - PIDMapT& particleIDMap); + template + std::vector + convertReconstructedParticles(const std::string& name, EVENT::LCCollection* LCCollection, RecoMapT& recoparticlesMap); /** * Convert a Vertex collection and return the resulting collection. @@ -286,17 +298,10 @@ namespace LCIO2EDM4hepConv { /** * Convert a Cluster collection and return the resulting collection. * Simultaneously populates the mapping from LCIO to EDM4hep objects. - * - * NOTE: Also populates a ParticleID collection, as those are persisted as - * part of the Cluster collection in LCIO. The name of this collection is - * _particleIDs - */ - template - std::vector convertClusters( - const std::string& name, - EVENT::LCCollection* LCCollection, - ClusterMapT& clusterMap, - PIDMapT& particleIDMap); + */ + template + std::unique_ptr + convertClusters(const std::string& name, EVENT::LCCollection* LCCollection, ClusterMapT& clusterMap); /** * Create an EventHeaderCollection and fills it with the Metadata. diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp index 68e6733b..9998c896 100644 --- a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp +++ b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp @@ -1,5 +1,8 @@ #include "k4EDM4hep2LcioConv/MappingUtils.h" +#include +#include + namespace LCIO2EDM4hepConv { template void convertObjectParameters(LCIOType* lcioobj, podio::Frame& event) @@ -114,15 +117,22 @@ namespace LCIO2EDM4hepConv { return dest; } - template - std::vector convertReconstructedParticles( - const std::string& name, - EVENT::LCCollection* LCCollection, - RecoMapT& recoparticlesMap, - PIDMapT& particleIDMap) + template + std::vector + convertReconstructedParticles(const std::string& name, EVENT::LCCollection* LCCollection, RecoMapT& recoparticlesMap) { auto dest = std::make_unique(); - auto particleIDs = std::make_unique(); + + // Set up a PIDHandler to split off the ParticlID objects stored in the + // reconstructed particles into separate collections. Each algorithm id / + // name get's a separate collection + auto pidHandler = UTIL::PIDHandler(LCCollection); + // TODO: parameter names + std::map> particleIDs; + for (const auto id : pidHandler.getAlgorithmIDs()) { + particleIDs.emplace(id, std::make_unique()); + } + for (unsigned i = 0, N = LCCollection->getNumberOfElements(); i < N; ++i) { auto* rval = static_cast(LCCollection->getElementAt(i)); auto lval = dest->create(); @@ -146,39 +156,26 @@ namespace LCIO2EDM4hepConv { } // Need to convert the particle IDs here, since they are part of the reco - // particle collection in LCIO + // particle collection in LCIO. for (const auto lcioPid : rval->getParticleIDs()) { auto pid = convertParticleID(lcioPid); - const auto [pidIt, pidInserted] = k4EDM4hep2LcioConv::detail::mapInsert( - lcioPid, pid, particleIDMap, k4EDM4hep2LcioConv::detail::InsertMode::Checked); - if (pidInserted) { - lval.addToParticleIDs(pid); - particleIDs->push_back(pid); + pid.setParticle(lval); + if (auto pidIt = particleIDs.find(pid.getAlgorithmType()); pidIt != particleIDs.end()) { + pidIt->second->push_back(pid); } else { - lval.addToParticleIDs(k4EDM4hep2LcioConv::detail::getMapped(pidIt)); + std::cerr << "ERROR: Found a PID object with an algorithm ID that is not known to the PIDHandler (id = " + << pid.getAlgorithmType() << ")" << std::endl; } } - - const auto lcioPidUsed = rval->getParticleIDUsed(); - if (lcioPidUsed == nullptr) { - continue; - } - if (const auto edm4hepPid = k4EDM4hep2LcioConv::detail::mapLookupTo(lcioPidUsed, particleIDMap)) { - lval.setParticleIDUsed(edm4hepPid.value()); - } - else { - auto pid = convertParticleID(lcioPidUsed); - particleIDs->push_back(pid); - k4EDM4hep2LcioConv::detail::mapInsert(lcioPidUsed, pid, particleIDMap); - lval.setParticleIDUsed(pid); - } } std::vector results; - results.reserve(2); + results.reserve(particleIDs.size() + 1); results.emplace_back(name, std::move(dest)); - results.emplace_back(name + "_particleIDs", std::move(particleIDs)); + for (auto& [id, coll] : particleIDs) { + results.emplace_back(getPIDCollName(name, pidHandler.getAlgorithmName(id)), std::move(coll)); + } return results; } @@ -472,14 +469,10 @@ namespace LCIO2EDM4hepConv { return dest; } - template - std::vector convertClusters( - const std::string& name, - EVENT::LCCollection* LCCollection, - ClusterMapT& clusterMap, - PIDMapT& particleIDMap) + template + std::unique_ptr + convertClusters(const std::string& name, EVENT::LCCollection* LCCollection, ClusterMapT& clusterMap) { - auto particleIDs = std::make_unique(); auto dest = std::make_unique(); for (unsigned i = 0, N = LCCollection->getNumberOfElements(); i < N; ++i) { @@ -503,26 +496,8 @@ namespace LCIO2EDM4hepConv { std::cerr << "EDM4hep element " << existingId << " did not get inserted. It belongs to the " << name << " collection" << std::endl; } - - // Need to convert the particle IDs here, since they are part of the cluster - // collection in LCIO - for (const auto lcioPid : rval->getParticleIDs()) { - auto pid = convertParticleID(lcioPid); - const auto [pidIt, pidInserted] = k4EDM4hep2LcioConv::detail::mapInsert(lcioPid, pid, particleIDMap); - if (pidInserted) { - lval.addToParticleIDs(pid); - particleIDs->push_back(pid); - } - else { - lval.addToParticleIDs(k4EDM4hep2LcioConv::detail::getMapped(pidIt)); - } - } } - std::vector results; - results.reserve(2); - results.emplace_back(name, std::move(dest)); - results.emplace_back(name + "_particleIDs", std::move(particleIDs)); - return results; + return dest; } template @@ -535,7 +510,7 @@ namespace LCIO2EDM4hepConv { retColls.emplace_back(name, convertMCParticles(name, LCCollection, typeMapping.mcParticles)); } else if (type == "ReconstructedParticle") { - return convertReconstructedParticles(name, LCCollection, typeMapping.recoParticles, typeMapping.particleIDs); + return convertReconstructedParticles(name, LCCollection, typeMapping.recoParticles); } else if (type == "Vertex") { retColls.emplace_back(name, convertVertices(name, LCCollection, typeMapping.vertices)); @@ -544,7 +519,7 @@ namespace LCIO2EDM4hepConv { retColls.emplace_back(name, convertTracks(name, LCCollection, typeMapping.tracks)); } else if (type == "Cluster") { - return convertClusters(name, LCCollection, typeMapping.clusters, typeMapping.particleIDs); + retColls.emplace_back(name, convertClusters(name, LCCollection, typeMapping.clusters)); } else if (type == "SimCalorimeterHit") { retColls.emplace_back(name, convertSimCalorimeterHits(name, LCCollection, typeMapping.simCaloHits)); diff --git a/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp b/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp index a32bedb3..40a3c4b2 100644 --- a/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp +++ b/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp @@ -1,6 +1,13 @@ #include "k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h" + #include "edm4hep/Constants.h" -#include "EVENT/MCParticle.h" +#include "edm4hep/utils/ParticleIDUtils.h" + +#include "UTIL/PIDHandler.h" +#include + +#include +#include namespace EDM4hep2LCIOConv { @@ -29,11 +36,51 @@ namespace EDM4hep2LCIOConv { return std::find(coll->begin(), coll->end(), collection_name) != coll->end(); } + void sortParticleIDs(std::vector& pidCollections) + { + std::sort(pidCollections.begin(), pidCollections.end(), [](const auto& pid1, const auto& pid2) { + static const auto defaultPidMeta = edm4hep::utils::ParticleIDMeta {"", std::numeric_limits::max(), {}}; + return pid1.metadata.value_or(defaultPidMeta).algoType < pid2.metadata.value_or(defaultPidMeta).algoType; + }); + } + + std::optional attachParticleIDMetaData( + IMPL::LCEventImpl* lcEvent, + const podio::Frame& edmEvent, + const ParticleIDConvData& pidCollMetaInfo) + { + const auto& [name, coll, pidMetaInfo] = pidCollMetaInfo; + const auto recoName = edmEvent.getName((*coll)[0].getParticle().id().collectionID); + // If we can't get the reconstructed particle collection name there is not + // much we can do + if (!recoName.has_value()) { + return std::nullopt; + } + // If we can't get meta data information there is not much we can do either + if (!pidMetaInfo.has_value()) { + return std::nullopt; + } + if (pidMetaInfo.has_value() && !recoName.has_value()) { + return pidMetaInfo->algoType; + } + + UTIL::PIDHandler pidHandler(lcEvent->getCollection(recoName.value())); + return pidHandler.addAlgorithm(pidMetaInfo->algoName, pidMetaInfo->paramNames); + } + std::unique_ptr convertEvent(const podio::Frame& edmEvent, const podio::Frame& metadata) { auto lcioEvent = std::make_unique(); auto objectMappings = CollectionsPairVectors {}; + // We have to convert these after all other (specifically + // ReconstructedParticle) collections have been converted. Otherwise we will + // not be able to set all the metadata for the PIDHandler (LCIO) to work + // properly. Here we store the name, the collection as well as potentially + // available meta information that we obtain when we first come across a + // collection + std::vector pidCollections {}; + const auto& collections = edmEvent.getAvailableCollections(); for (const auto& name : collections) { const auto edmCollection = edmEvent.get(name); @@ -92,11 +139,11 @@ namespace EDM4hep2LCIOConv { else if (auto coll = dynamic_cast(edmCollection)) { convertEventHeader(coll, lcioEvent.get()); } - else if ( - dynamic_cast(edmCollection) || - dynamic_cast(edmCollection)) { - // "converted" as part of FillMissingCollectoins at the end or as part - // of the reconstructed particle + else if (auto coll = dynamic_cast(edmCollection)) { + pidCollections.emplace_back(name, coll, edm4hep::utils::PIDHandler::getAlgoInfo(metadata, name)); + } + else if (dynamic_cast(edmCollection)) { + // "converted" during relation resolving later continue; } else { @@ -110,6 +157,14 @@ namespace EDM4hep2LCIOConv { } } + sortParticleIDs(pidCollections); + for (const auto& pidCollMeta : pidCollections) { + // Use -1 as a somewhat easy to identify value of missing reco collections + // or pid metadata + const auto algoId = attachParticleIDMetaData(lcioEvent.get(), edmEvent, pidCollMeta).value_or(-1); + convertParticleIDs(pidCollMeta.coll, objectMappings.particleIDs, algoId); + } + resolveRelations(objectMappings); return lcioEvent; diff --git a/k4EDM4hep2LcioConv/src/k4Lcio2EDM4hepConv.cpp b/k4EDM4hep2LcioConv/src/k4Lcio2EDM4hepConv.cpp index a87ad03d..d49ede84 100644 --- a/k4EDM4hep2LcioConv/src/k4Lcio2EDM4hepConv.cpp +++ b/k4EDM4hep2LcioConv/src/k4Lcio2EDM4hepConv.cpp @@ -1,5 +1,7 @@ #include "k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.h" +#include + #include namespace LCIO2EDM4hepConv { @@ -71,6 +73,17 @@ namespace LCIO2EDM4hepConv { return headerColl; } + std::vector getPIDMetaInfo(const EVENT::LCCollection* recoColl) + { + std::vector pidInfos {}; + const auto pidHandler = UTIL::PIDHandler(recoColl); + for (const auto id : pidHandler.getAlgorithmIDs()) { + pidInfos.emplace_back(pidHandler.getAlgorithmName(id), id, pidHandler.getParameterNames(id)); + } + + return pidInfos; + } + podio::Frame convertEvent(EVENT::LCEvent* evt, const std::vector& collsToConvert) { auto typeMapping = LcioEdmTypeMapping {}; diff --git a/standalone/lcio2edm4hep.cpp b/standalone/lcio2edm4hep.cpp index 496bde3d..bff47177 100644 --- a/standalone/lcio2edm4hep.cpp +++ b/standalone/lcio2edm4hep.cpp @@ -4,6 +4,8 @@ #include #include +#include + #include "podio/podioVersion.h" #if PODIO_BUILD_VERSION >= PODIO_VERSION(0, 99, 0) #include "podio/ROOTWriter.h" @@ -166,6 +168,8 @@ int main(int argc, char* argv[]) podio::ROOTWriter writer(args.outputFile); + podio::Frame metadata {}; + for (int j = 0; j < lcreader->getNumberOfRuns(); ++j) { if (j % 1 == 0) { std::cout << "processing RunHeader: " << j << std::endl; @@ -187,9 +191,26 @@ int main(int argc, char* argv[]) colPatcher.patchCollections(evt); } const auto edmEvent = LCIO2EDM4hepConv::convertEvent(evt, collsToConvert); + + // For the first event we also convert some meta information for the + // ParticleID handling + if (i == 0) { + for (const auto& name : *evt->getCollectionNames()) { + auto coll = evt->getCollection(name); + if (coll->getTypeName() == "ReconstructedParticle") { + for (const auto& pidInfo : LCIO2EDM4hepConv::getPIDMetaInfo(coll)) { + edm4hep::utils::PIDHandler::setAlgoInfo( + metadata, LCIO2EDM4hepConv::getPIDCollName(name, pidInfo.algoName), pidInfo); + } + } + } + } + writer.writeFrame(edmEvent, "events"); } + writer.writeFrame(metadata, podio::Category::Metadata); + writer.finish(); return 0; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ae868cc..b275acdf 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,9 @@ add_library(edmCompare SHARED src/CompareEDM4hepLCIO.cc src/ObjectMapping.cc src/CompareEDM4hepEDM4hep.cc) target_link_libraries(edmCompare PUBLIC EDM4HEP::edm4hep LCIO::lcio) -target_include_directories(edmCompare PUBLIC ${LCIO_INCLUDE_DIRS}) +target_include_directories(edmCompare + PUBLIC + $ +) add_library(TestUtils SHARED src/EDM4hep2LCIOUtilities.cc) target_link_libraries(TestUtils PUBLIC EDM4HEP::edm4hep LCIO::lcio) diff --git a/tests/edm4hep_roundtrip.cpp b/tests/edm4hep_roundtrip.cpp index 0438a6b4..2ed4f24f 100644 --- a/tests/edm4hep_roundtrip.cpp +++ b/tests/edm4hep_roundtrip.cpp @@ -16,8 +16,8 @@ int main() { - const auto origEvent = createExampleEvent(); - const auto lcioEvent = EDM4hep2LCIOConv::convertEvent(origEvent); + const auto& [origEvent, metadata] = createExampleEvent(); + const auto lcioEvent = EDM4hep2LCIOConv::convertEvent(origEvent, metadata); const auto roundtripEvent = LCIO2EDM4hepConv::convertEvent(lcioEvent.get()); ASSERT_SAME_OR_ABORT(edm4hep::CalorimeterHitCollection, "caloHits"); @@ -28,6 +28,9 @@ int main() ASSERT_SAME_OR_ABORT(edm4hep::TrackerHitPlaneCollection, "trackerHitPlanes"); ASSERT_SAME_OR_ABORT(edm4hep::ClusterCollection, "clusters"); ASSERT_SAME_OR_ABORT(edm4hep::ReconstructedParticleCollection, "recos"); + ASSERT_SAME_OR_ABORT(edm4hep::ParticleIDCollection, "ParticleID_coll_1"); + ASSERT_SAME_OR_ABORT(edm4hep::ParticleIDCollection, "ParticleID_coll_2"); + ASSERT_SAME_OR_ABORT(edm4hep::ParticleIDCollection, "ParticleID_coll_3"); return 0; } diff --git a/tests/edm4hep_to_lcio.cpp b/tests/edm4hep_to_lcio.cpp index fb072dd2..22cc9f91 100644 --- a/tests/edm4hep_to_lcio.cpp +++ b/tests/edm4hep_to_lcio.cpp @@ -11,8 +11,8 @@ int main() { - const auto edmEvent = createExampleEvent(); - const auto lcioEvent = EDM4hep2LCIOConv::convertEvent(edmEvent); + const auto& [edmEvent, metadata] = createExampleEvent(); + const auto lcioEvent = EDM4hep2LCIOConv::convertEvent(edmEvent, metadata); if (!compareEventHeader(lcioEvent.get(), &edmEvent)) { return 1; diff --git a/tests/src/CompareEDM4hepEDM4hep.cc b/tests/src/CompareEDM4hepEDM4hep.cc index ae2ddea7..05b54256 100644 --- a/tests/src/CompareEDM4hepEDM4hep.cc +++ b/tests/src/CompareEDM4hepEDM4hep.cc @@ -8,6 +8,7 @@ #include "edm4hep/TrackerHitPlaneCollection.h" #include "edm4hep/ClusterCollection.h" #include "edm4hep/ReconstructedParticleCollection.h" +#include "edm4hep/ParticleIDCollection.h" #include #include @@ -321,17 +322,32 @@ bool compare( relParticles[iP].getObjectID(), "related particle " << iP << " in reco particle " << i); } + } - const auto origRelPIDs = origReco.getParticleIDs(); - const auto relPIDs = reco.getParticleIDs(); - REQUIRE_SAME(origRelPIDs.size(), relPIDs.size(), "number of related ParticleIDs in reco particle" << i); - for (size_t iP = 0; iP < relPIDs.size(); ++iP) { - REQUIRE_SAME( - origRelPIDs[iP].getObjectID(), - relPIDs[iP].getObjectID(), - "related ParticleID " << iP << " in reco particle " << i); + return true; +} + +bool compare(const edm4hep::ParticleIDCollection& origColl, const edm4hep::ParticleIDCollection& roundtripColl) +{ + REQUIRE_SAME(origColl.size(), roundtripColl.size(), "collection sizes"); + for (size_t i = 0; i < origColl.size(); ++i) { + const auto origPid = origColl[i]; + const auto pid = roundtripColl[i]; + + // This might not be preserved in roundtripping + REQUIRE_SAME(origPid.getAlgorithmType(), pid.getAlgorithmType(), "algorithm type in ParticleID " << i); + + REQUIRE_SAME(origPid.getType(), pid.getType(), "type in ParticleID " << i); + + const auto origParams = origPid.getParameters(); + const auto params = pid.getParameters(); + REQUIRE_SAME(origParams.size(), params.size(), "parameter sizes in ParticleID " << i); + for (size_t iP = 0; iP < params.size(); ++iP) { + REQUIRE_SAME(origParams[iP], params[iP], "parameter " << iP << " in ParticleID " << i); } - } + REQUIRE_SAME( + origPid.getParticle().getObjectID(), pid.getParticle().getObjectID(), "related particle in ParticleID " << i); + } return true; } diff --git a/tests/src/CompareEDM4hepEDM4hep.h b/tests/src/CompareEDM4hepEDM4hep.h index 20d6860c..7d2f1881 100644 --- a/tests/src/CompareEDM4hepEDM4hep.h +++ b/tests/src/CompareEDM4hepEDM4hep.h @@ -25,4 +25,6 @@ bool compare( const edm4hep::ReconstructedParticleCollection& origColl, const edm4hep::ReconstructedParticleCollection& roundtripColl); +bool compare(const edm4hep::ParticleIDCollection& origColl, const edm4hep::ParticleIDCollection& roundtripColl); + #endif // K4EDM4HEP2LCIOCONV_TEST_COMPAREEDM4HEPEDM4HEP_H diff --git a/tests/src/CompareEDM4hepLCIO.cc b/tests/src/CompareEDM4hepLCIO.cc index 5a6d4e14..08160308 100644 --- a/tests/src/CompareEDM4hepLCIO.cc +++ b/tests/src/CompareEDM4hepLCIO.cc @@ -18,8 +18,6 @@ * the single object compare. All comparison functions return true if all * comparisons succeeded or false in case any comparison failed returning as * early as possible. - * - * TODO: Also compare relations */ /// Convert the two 32 bit cellIDs into one 64 bit value @@ -193,34 +191,26 @@ bool compare( ASSERT_COMPARE_RELATION( lcioElem, edm4hepElem, getStartVertex, objectMaps.vertices, "startVertex in ReconstructedParticle"); + // ParticleIDs need special treatment because they live in different + // collections in EDM4hep. Here we make sure that all ParticleIDs have been + // converted and mapped correctly by checking the ParticleID content and + // making sure that it points back to the converted reco particle const auto& lcioPIDs = lcioElem->getParticleIDs(); - const auto edmPIDs = edm4hepElem.getParticleIDs(); - ASSERT_COMPARE_VALS(lcioPIDs.size(), edmPIDs.size(), "particleIDs with different sizes in ReconstructedParticle"); - for (size_t i = 0; i < lcioPIDs.size(); ++i) { - if (!compare(lcioPIDs[i], edmPIDs[i])) { - std::cerr << "particle ID " << i << " differs in ReconstructedParticle (LCIO: " << lcioPIDs[i] - << ", EDM4hep: " << edmPIDs[i] << ")" << std::endl; - return false; - } - } - - const auto lcioPIDUsed = lcioElem->getParticleIDUsed(); - const auto edmPIDUsed = edm4hepElem.getParticleIDUsed(); - if (lcioPIDUsed == nullptr) { - if (edmPIDUsed.isAvailable()) { - std::cerr << "particleIDUsed is not available in LCIO, but points to " << edmPIDUsed.getObjectID() - << " in EDM4hep for ReconstructedParticle" << std::endl; - return false; + if (auto it = objectMaps.particleIDs.find(lcioPIDs[i]); it != objectMaps.particleIDs.end()) { + const auto& [lcioPid, edm4hepPid] = *it; + if (!compare(lcioPid, edm4hepPid) || edm4hepPid.getParticle() != edm4hepElem) { + std::cerr << "particle ID " << i << " is not mapped to the correct EDM4hep particle ID (LCIO: " << lcioPid + << ", EDM4hep: " << edm4hepPid << ")" << std::endl; + return false; + } } - } - else { - if (!compare(lcioPIDUsed, edmPIDUsed)) { - std::cerr << "particleIDUsed differs in ReconstructedParticle (LCIO: " << lcioPIDUsed - << ", EDM4hep: " << edmPIDUsed << ")" << std::endl; + else { + std::cerr << "Cannot find a converted ParticleID object for particle ID " << i << std::endl; return false; } } + return true; } diff --git a/tests/src/CompareEDM4hepLCIO.h b/tests/src/CompareEDM4hepLCIO.h index 977496d1..22f05c59 100644 --- a/tests/src/CompareEDM4hepLCIO.h +++ b/tests/src/CompareEDM4hepLCIO.h @@ -34,6 +34,8 @@ namespace edm4hep { #endif #include "edm4hep/TrackerHitPlaneCollection.h" #include "edm4hep/VertexCollection.h" +#include "edm4hep/utils/ParticleIDUtils.h" + #include "podio/Frame.h" #include @@ -53,6 +55,7 @@ namespace edm4hep { #include #include #include +#include #include bool compare( diff --git a/tests/src/EDM4hep2LCIOUtilities.cc b/tests/src/EDM4hep2LCIOUtilities.cc index ca6507d0..ff43a808 100644 --- a/tests/src/EDM4hep2LCIOUtilities.cc +++ b/tests/src/EDM4hep2LCIOUtilities.cc @@ -21,6 +21,7 @@ namespace edm4hep { #include "edm4hep/ClusterCollection.h" #include "edm4hep/ReconstructedParticleCollection.h" #include +#include #include "podio/Frame.h" @@ -328,17 +329,15 @@ edm4hep::ClusterCollection createClusters( return clusterColl; } -std::tuple createRecoParticles( +edm4hep::ReconstructedParticleCollection createRecoParticles( const int nRecos, const edm4hep::TrackCollection& tracks, const std::vector& trackIdcs, const edm4hep::ClusterCollection& clusters, const std::vector& clusterIdcs, - const std::vector& recIdcs, - const std::vector& pidAlgTypes) + const std::vector& recIdcs) { auto recoColl = edm4hep::ReconstructedParticleCollection {}; - auto pidColl = edm4hep::ParticleIDCollection {}; for (int i = 0; i < nRecos; ++i) { auto reco = recoColl.create(); reco.setCharge(1.23f); @@ -355,18 +354,35 @@ std::tuple createParticleIDs( + const std::vector>& recoIdcs, + const edm4hep::ReconstructedParticleCollection& recoParticles) +{ + std::vector collections; + collections.reserve(recoIdcs.size()); + + float param = 0; + for (const auto& idcs : recoIdcs) { + auto& coll = collections.emplace_back(); + for (const auto idx : idcs) { + auto pid = coll.create(); + pid.setType(idx); + pid.setParticle(recoParticles[idx]); + pid.addToParameters(param++); + } } - return {std::move(recoColl), std::move(pidColl)}; + return collections; } -podio::Frame createExampleEvent() +std::tuple createExampleEvent() { - podio::Frame event; + auto retTuple = std::make_tuple(podio::Frame {}, podio::Frame {}); + + auto& [event, metadata] = retTuple; event.put(createEventHeader(), "EventHeader"); const auto& mcParticles = @@ -401,18 +417,29 @@ podio::Frame createExampleEvent() event.put(std::move(tmpSimCaloHits), "simCaloHits"); event.put(std::move(tmpCaloHitConts), "caloHitContributions"); - auto [recColl, pidColl] = createRecoParticles( - test_config::nRecoParticles, - tracks, - test_config::recoTrackIdcs, - clusters, - test_config::recoClusterIdcs, - test_config::recoRecoIdcs, - test_config::recoPIDTypes); - - event.put(std::move(recColl), "recos"); - // Make sure to use the same name as is generated for the LCIO to EDM4hep conversion - event.put(std::move(pidColl), "recos_particleIDs"); + const auto& recoColl = event.put( + createRecoParticles( + test_config::nRecoParticles, + tracks, + test_config::recoTrackIdcs, + clusters, + test_config::recoClusterIdcs, + test_config::recoRecoIdcs), + "recos"); + + // Start at 0 here because that is also where the PIDHandler in LCIO starts + int algoId = 0; + for (auto& pidColl : createParticleIDs(test_config::pidRecoIdcs, recoColl)) { + // Make sure to use the same name as is generated for the LCIO to EDM4hep + // conversion + const auto pidCollName = "recos_PID_pidAlgo_" + std::to_string(algoId); + edm4hep::utils::PIDHandler::setAlgoInfo( + metadata, pidColl, pidCollName, {"pidAlgo_" + std::to_string(algoId), algoId, {"param"}}); + + event.put(std::move(pidColl), pidCollName); + + algoId++; + } - return event; + return retTuple; } diff --git a/tests/src/EDM4hep2LCIOUtilities.h b/tests/src/EDM4hep2LCIOUtilities.h index 210256af..94ae10d5 100644 --- a/tests/src/EDM4hep2LCIOUtilities.h +++ b/tests/src/EDM4hep2LCIOUtilities.h @@ -108,11 +108,8 @@ namespace test_config { /// particle const static std::vector recoRecoIdcs = {{0, 3}, {2, 2}, {5, 1}, {5, 0}, {1, 2}, {2, 4}}; - /// The ParticleIDs algorithmType that will be create for some reco particles. - /// The first index is the reco paritcle to which a ParticleID with algorithm - /// type of the second index will be attached - const static std::vector recoPIDTypes = {{0, 1}, {0, 2}, {1, 1}, {2, 1}, {2, 2}, {3, 1}, {4, 2}}; - + /// The number of entries for the generated ParticleID collections + const static std::vector> pidRecoIdcs = {{1, 3, 4}, {2, 3}, {0, 1, 2, 3, 4, 5}}; } // namespace test_config /** @@ -176,17 +173,21 @@ edm4hep::ClusterCollection createClusters( const std::vector& clusterHitIdcs, const std::vector& clusterClusterIdcs); -std::tuple createRecoParticles( +edm4hep::ReconstructedParticleCollection createRecoParticles( const int nRecos, const edm4hep::TrackCollection& tracks, const std::vector& trackIdcs, const edm4hep::ClusterCollection& clusters, const std::vector& clusterIdcs, - const std::vector& recIdcs, - const std::vector& pidAlgTypes); + const std::vector& recIdcs); + +std::vector createParticleIDs( + const std::vector>& recoIdcs, + const edm4hep::ReconstructedParticleCollection& recoParticles); /** - * Create an example event that can be used to test the converter. + * Create an example event that can be used to test the converter. Also populate + * a metadata frame that can be used in the conversion. * * Content: * @@ -207,6 +208,6 @@ std::tuple createExampleEvent(); #endif // K4EDM4HEP2LCIOCONV_TEST_COMPAREEDM4HEPLCIO_H diff --git a/tests/src/ObjectMapping.cc b/tests/src/ObjectMapping.cc index a6ac8bb8..3a7b81b4 100644 --- a/tests/src/ObjectMapping.cc +++ b/tests/src/ObjectMapping.cc @@ -15,7 +15,9 @@ #include "EVENT/MCParticle.h" #include "EVENT/LCEvent.h" #include "EVENT/LCCollection.h" +#include "EVENT/ParticleID.h" #include "UTIL/LCIterator.h" +#include "UTIL/PIDHandler.h" #include "edm4hep/TrackCollection.h" #if __has_include("edm4hep/TrackerHit3DCollection.h") @@ -36,9 +38,16 @@ namespace edm4hep { #include "edm4hep/ReconstructedParticleCollection.h" #include "edm4hep/RawTimeSeriesCollection.h" #include "edm4hep/VertexCollection.h" +#include "edm4hep/ParticleIDCollection.h" #include "podio/Frame.h" +#include "k4EDM4hep2LcioConv/MappingUtils.h" + +#include +#include +#include + template void fillMap(MapT& map, EVENT::LCCollection* lcioColl, const EDM4hepT& edmColl) { @@ -58,10 +67,46 @@ void fillMap(MapT& map, EVENT::LCCollection* lcioColl, const EDM4hepT& edmColl) fillMap(mapName, lcioColl, edm4hepColl); \ } +void fillRecoPIDMaps( + ObjectMappings::Map& recoMap, + std::unordered_map& pidMap, + const std::string& recoName, + EVENT::LCEvent* lcEvt, + const podio::Frame& edmEvt) +{ + UTIL::LCIterator lcioIt(lcEvt, recoName); + UTIL::PIDHandler pidHandler(lcioIt()); + const auto& edmColl = edmEvt.get(recoName); + + // Simply need to assume here that per algorithm the pid "collections" run in + // parallel. + std::unordered_map algoIdcs {}; + for (const auto edmElem : edmColl) { + const auto* lcioElem = lcioIt.next(); + recoMap.emplace(lcioElem, edmElem.getObjectID()); + + for (const auto* lcioPid : lcioElem->getParticleIDs()) { + const auto& algoName = pidHandler.getAlgorithmName(lcioPid->getAlgorithmType()); + auto idx = algoIdcs[algoName]++; + // No need to do any fancy caching, because the Frame does that already + const auto& pidColl = edmEvt.get(recoName + "_PID_" + algoName); + pidMap.emplace(lcioPid, pidColl[idx]); + } + } +} + +std::string getRecoName(const std::string& pidName) +{ + const auto pos = pidName.find("_PID_"); + if (pos != std::string::npos) { + return pidName.substr(0, pos); + } + return ""; +} + ObjectMappings ObjectMappings::fromEvent(EVENT::LCEvent* lcEvt, const podio::Frame& edmEvt) { ObjectMappings mapping {}; - for (const auto& name : *(lcEvt->getCollectionNames())) { // We only use non subset collections here, because we want the "real" // objects @@ -81,7 +126,6 @@ ObjectMappings ObjectMappings::fromEvent(EVENT::LCEvent* lcEvt, const podio::Fra FILL_MAP(RawCalorimeterHit, mapping.rawCaloHits); FILL_MAP(SimCalorimeterHit, mapping.simCaloHits); FILL_MAP(MCParticle, mapping.mcParticles); - FILL_MAP(ReconstructedParticle, mapping.recoParticles); FILL_MAP(Vertex, mapping.vertices); // Need special treatment for TPCHit type mismatch if (type == "TPCHit") { @@ -92,6 +136,11 @@ ObjectMappings ObjectMappings::fromEvent(EVENT::LCEvent* lcEvt, const podio::Fra auto& edm4hepColl = edmEvt.get(name); fillMap(mapping.trackerHits, lcioColl, edm4hepColl); } + // Special treatment for reco particles and ParticleIDs because of + // conceptual differences + if (type == "ReconstructedParticle") { + fillRecoPIDMaps(mapping.recoParticles, mapping.particleIDs, name, lcEvt, edmEvt); + } } return mapping; diff --git a/tests/src/ObjectMapping.h b/tests/src/ObjectMapping.h index 00bd04b7..9e472f17 100644 --- a/tests/src/ObjectMapping.h +++ b/tests/src/ObjectMapping.h @@ -18,12 +18,17 @@ namespace EVENT { class ReconstructedParticle; class MCParticle; class LCEvent; + class ParticleID; } // namespace EVENT namespace podio { class Frame; } // namespace podio +namespace edm4hep { + class ParticleID; +} + struct ObjectMappings { template using Map = std::unordered_map; @@ -41,6 +46,8 @@ struct ObjectMappings { Map tpcHits {}; Map vertices {}; + std::unordered_map particleIDs; + static ObjectMappings fromEvent(EVENT::LCEvent* lcEvt, const podio::Frame& edmEvt); };