EDM4hep/src/EDM4hepUtil.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp"
0010
0011 #include "Acts/Definitions/Common.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
0016 #include "ActsExamples/EventData/Index.hpp"
0017 #include "ActsExamples/EventData/Measurement.hpp"
0018 #include "ActsExamples/Validation/TrackClassification.hpp"
0019 #include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp"
0020 #include "ActsPlugins/EDM4hep/EDM4hepUtil.hpp"
0021 #include <ActsPodioEdm/TrackerHitLocal.h>
0022
0023 #include <cstdint>
0024
0025 #include "edm4hep/TrackState.h"
0026
0027 using namespace Acts::UnitLiterals;
0028
0029 namespace ActsExamples {
0030
0031 SimParticle EDM4hepUtil::readParticle(const edm4hep::MCParticle& from,
0032                                       const MapParticleIdFrom& particleMapper) {
0033   SimBarcode particleId = particleMapper(from);
0034
0035   SimParticle to(particleId, static_cast<Acts::PdgParticle>(from.getPDG()),
0036                  from.getCharge() * Acts::UnitConstants::e,
0037                  from.getMass() * Acts::UnitConstants::GeV);
0038
0039   // TODO do we have that in EDM4hep?
0040   // particle.setProcess(static_cast<ActsFatras::ProcessType>(data.process));
0041
0042   to.initialState().setPosition4(from.getVertex()[0] * Acts::UnitConstants::mm,
0043                                  from.getVertex()[1] * Acts::UnitConstants::mm,
0044                                  from.getVertex()[2] * Acts::UnitConstants::mm,
0045                                  from.getTime() * Acts::UnitConstants::ns);
0046
0047   // Only used for direction; normalization/units do not matter
0048   Acts::Vector3 momentum = {from.getMomentum()[0], from.getMomentum()[1],
0049                             from.getMomentum()[2]};
0050   to.initialState().setDirection(momentum.normalized());
0051
0052   to.initialState().setAbsoluteMomentum(momentum.norm() * 1_GeV);
0053
0054   return to;
0055 }
0056
0057 void EDM4hepUtil::writeParticle(const SimParticle& from,
0058                                 edm4hep::MutableMCParticle to) {
0059   // TODO what about particleId?
0060
0061   to.setPDG(from.pdg());
0062   to.setCharge(from.charge() / Acts::UnitConstants::e);
0063   to.setMass(from.mass() / Acts::UnitConstants::GeV);
0064   to.setVertex({from.position().x(), from.position().y(), from.position().z()});
0065   to.setMomentum({static_cast<float>(from.fourMomentum().x()),
0066                   static_cast<float>(from.fourMomentum().y()),
0067                   static_cast<float>(from.fourMomentum().z())});
0068   to.setMomentumAtEndpoint(
0069       {static_cast<float>(from.finalState().fourMomentum().x()),
0070        static_cast<float>(from.finalState().fourMomentum().y()),
0071        static_cast<float>(from.finalState().fourMomentum().z())});
0072 }
0073
0074 ActsFatras::Hit EDM4hepUtil::readSimHit(const edm4hep::SimTrackerHit& from,
0075                                         const MapParticleIdFrom& particleMapper,
0076                                         const MapGeometryIdFrom& geometryMapper,
0077                                         std::uint32_t index) {
0078   auto particle = ActsPlugins::EDM4hepUtil::getParticle(from);
0079   SimBarcode particleId = particleMapper(particle);
0080
0081   const auto mass = particle.getMass() * 1_GeV;
0082   const Acts::Vector3 momentum{
0083       from.getMomentum().x * 1_GeV,
0084       from.getMomentum().y * 1_GeV,
0085       from.getMomentum().z * 1_GeV,
0086   };
0087   const auto energy = std::hypot(momentum.norm(), mass);
0088
0089   Acts::Vector4 pos4{
0090       from.getPosition().x * 1_mm,
0091       from.getPosition().y * 1_mm,
0092       from.getPosition().z * 1_mm,
0093       from.getTime() * 1_ns,
0094   };
0095
0096   Acts::Vector4 mom4{
0097       momentum.x(),
0098       momentum.y(),
0099       momentum.z(),
0100       energy,
0101   };
0102
0103   Acts::Vector4 delta4 = Acts::Vector4::Zero();
0104   delta4[Acts::eEnergy] = -from.getEDep() * Acts::UnitConstants::GeV;
0105
0106   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0107
0108   return ActsFatras::Hit(geometryId, particleId, pos4, mom4, mom4 + delta4,
0109                          index);
0110 }
0111
0112 void EDM4hepUtil::writeSimHit(const ActsFatras::Hit& from,
0113                               edm4hep::MutableSimTrackerHit to,
0114                               const MapParticleIdTo& particleMapper,
0115                               const MapGeometryIdTo& geometryMapper) {
0116   const Acts::Vector4& globalPos4 = from.fourPosition();
0117   const Acts::Vector4& momentum4Before = from.momentum4Before();
0118   const auto delta4 = from.momentum4After() - momentum4Before;
0119
0120   if (particleMapper) {
0121     ActsPlugins::EDM4hepUtil::setParticle(to,
0122                                           particleMapper(from.particleId()));
0123   }
0124
0125   if (geometryMapper) {
0126     // TODO what about the digitization?
0127     to.setCellID(geometryMapper(from.geometryId()));
0128   }
0129
0130   to.setTime(globalPos4[Acts::eTime] / Acts::UnitConstants::ns);
0131
0132   to.setPosition({
0133       globalPos4[Acts::ePos0] / Acts::UnitConstants::mm,
0134       globalPos4[Acts::ePos1] / Acts::UnitConstants::mm,
0135       globalPos4[Acts::ePos2] / Acts::UnitConstants::mm,
0136   });
0137
0138   to.setMomentum({
0139       static_cast<float>(momentum4Before[Acts::eMom0] /
0140                          Acts::UnitConstants::GeV),
0141       static_cast<float>(momentum4Before[Acts::eMom1] /
0142                          Acts::UnitConstants::GeV),
0143       static_cast<float>(momentum4Before[Acts::eMom2] /
0144                          Acts::UnitConstants::GeV),
0145   });
0146
0147   to.setEDep(-delta4[Acts::eEnergy] / Acts::UnitConstants::GeV);
0148 }
0149
0150 VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement(
0151     MeasurementContainer& container, const edm4hep::TrackerHitPlane& from,
0152     const edm4hep::TrackerHit3DCollection* /*fromClusters*/,
0153     Cluster* /*toCluster*/, const MapGeometryIdFrom& geometryMapper) {
0154   // no need for digitization as we only want to identify the sensor
0155   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0156
0157   auto pos = from.getPosition();
0158   auto cov = from.getCovMatrix();
0159
0160   DigitizedParameters dParameters;
0161
0162   dParameters.indices.push_back(Acts::eBoundLoc0);
0163   dParameters.values.push_back(pos.x);
0164   dParameters.variances.push_back(cov[0]);
0165
0166   // TODO cut this out for 1D
0167   dParameters.indices.push_back(Acts::eBoundLoc1);
0168   dParameters.values.push_back(pos.y);
0169   dParameters.variances.push_back(cov[2]);
0170
0171   dParameters.indices.push_back(Acts::eBoundTime);
0172   dParameters.values.push_back(pos.z);
0173   dParameters.variances.push_back(cov[5]);
0174
0175   auto to = createMeasurement(container, geometryId, dParameters);
0176
0177   // @TODO: Figure out if cell information is accessible
0178
0179   return to;
0180 }
0181
0182 void EDM4hepUtil::writeMeasurement(
0183     const Acts::GeometryContext& gctx,
0184     const ConstVariableBoundMeasurementProxy& from,
0185     ActsPodioEdm::MutableTrackerHitLocal& to, const Acts::Surface& surface) {
0186   long dim = from.size();
0187
0188   const auto* placement = surface.surfacePlacement();
0189   if (placement == nullptr) {
0190     throw std::runtime_error("Surface placement not found");
0191   }
0192   const auto* dd4hepDetectorElement =
0193       dynamic_cast<const ActsPlugins::DD4hepDetectorElement*>(placement);
0194   if (dd4hepDetectorElement == nullptr) {
0195     throw std::runtime_error(
0196         "Surface placement is not a DD4hepDetectorElement");
0197   }
0198
0199   std::uint64_t cellId = dd4hepDetectorElement->sourceElement().volumeID();
0200
0201   ActsPlugins::EDM4hepUtil::writeMeasurement(
0202       gctx, {from.parameters().data(), dim},
0203       {from.covariance().data(), dim, dim}, from.subspaceHelper().indices(),
0204       cellId, surface, to);
0205 }
0206
0207 VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement(
0208     MeasurementContainer& container, const ActsPodioEdm::TrackerHitLocal& from,
0209     const MapGeometryIdFrom& geometryMapper) {
0210   auto data = ActsPlugins::EDM4hepUtil::readMeasurement(from);
0211   Acts::GeometryIdentifier geometryId = geometryMapper(data.cellId);
0212
0213   const auto dim = data.indices.size();
0214   Acts::DynamicVector parameters(dim);
0215   Acts::DynamicMatrix covariance(dim, dim);
0216   for (std::size_t i = 0; i < dim; ++i) {
0217     parameters(i) = data.parameters(data.indices[i]);
0218     for (std::size_t j = 0; j < dim; ++j) {
0219       covariance(i, j) = data.covariance(data.indices[i], data.indices[j]);
0220     }
0221   }
0222
0223   return container.emplaceMeasurement(static_cast<std::uint8_t>(dim),
0224                                       geometryId, data.indices, parameters,
0225                                       covariance);
0226 }
0227
0228 void EDM4hepUtil::writeTrajectory(
0229     const Acts::GeometryContext& gctx, double Bz, const Trajectories& from,
0230     edm4hep::MutableTrack to, std::size_t fromIndex,
0231     const Acts::ParticleHypothesis& particleHypothesis,
0232     const MeasurementParticlesMap& measurementParticlesMap) {
0233   const auto& multiTrajectory = from.multiTrajectory();
0234   auto trajectoryState =
0235       Acts::MultiTrajectoryHelpers::trajectoryState(multiTrajectory, fromIndex);
0236
0237   std::vector<ParticleHitCount> particleHitCount;
0238   identifyContributingParticles(measurementParticlesMap, from, fromIndex,
0239                                 particleHitCount);
0240   // TODO use particles
0241
0242   // TODO write track params
0243   // auto trackParameters = from.trackParameters(fromIndex);
0244
0245   to.setChi2(trajectoryState.chi2Sum / trajectoryState.NDF);
0246   to.setNdf(trajectoryState.NDF);
0247
0248   multiTrajectory.visitBackwards(fromIndex, [&](const auto& state) {
0249     // we only fill the track states with non-outlier measurement
0250     auto typeFlags = state.typeFlags();
0251     if (!typeFlags.isMeasurement()) {
0252       return true;
0253     }
0254
0255     edm4hep::TrackState trackState;
0256
0257     Acts::BoundTrackParameters parObj{state.referenceSurface().getSharedPtr(),
0258                                       state.parameters(), state.covariance(),
0259                                       particleHypothesis};
0260
0261     // Convert to LCIO track parametrization expected by EDM4hep
0262     // This will create an ad-hoc perigee surface if the input parameters are
0263     // not bound on a perigee surface already
0264     ActsPlugins::EDM4hepUtil::detail::Parameters converted =
0265         ActsPlugins::EDM4hepUtil::detail::convertTrackParametersToEdm4hep(
0266             gctx, Bz, parObj);
0267
0268     trackState.D0 = converted.values[0];
0269     trackState.Z0 = converted.values[1];
0270     trackState.phi = converted.values[2];
0271     trackState.tanLambda = converted.values[3];
0272     trackState.omega = converted.values[4];
0273     trackState.time = converted.values[5];
0274
0275     // Converted parameters are relative to an ad-hoc perigee surface created
0276     // at the hit location
0277     auto center = converted.surface->center(gctx);
0278     trackState.referencePoint.x = center.x();
0279     trackState.referencePoint.y = center.y();
0280     trackState.referencePoint.z = center.z();
0281
0282     if (converted.covariance) {
0283       auto c = [&](std::size_t row, std::size_t col) {
0284         return static_cast<float>(converted.covariance.value()(row, col));
0285       };
0286
0287       // clang-format off
0288       trackState.covMatrix = edm4hep::CovMatrix6f{
0289         c(0, 0),
0290         c(1, 0), c(1, 1),
0291         c(2, 0), c(2, 1), c(2, 2),
0292         c(3, 0), c(3, 1), c(3, 2), c(3, 3),
0293         c(4, 0), c(4, 1), c(4, 2), c(4, 3), c(4, 4),
0294         c(5, 0), c(5, 1), c(5, 2), c(5, 3), c(5, 4), c(5, 5)
0295       };
0296       // clang-format on
0297     }
0298
0299     to.addToTrackStates(trackState);
0300
0301     return true;
0302   });
0303 }
0304
0305 }  // namespace ActsExamples