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/EDM4hep/EDM4hepUtil.hpp"
0020
0021 #include "edm4hep/TrackState.h"
0022
0023 using namespace Acts::UnitLiterals;
0024
0025 namespace ActsExamples {
0026
0027 SimParticle EDM4hepUtil::readParticle(const edm4hep::MCParticle& from,
0028                                       const MapParticleIdFrom& particleMapper) {
0029   ActsFatras::Barcode particleId = particleMapper(from);
0030
0031   SimParticle to(particleId, static_cast<Acts::PdgParticle>(from.getPDG()),
0032                  from.getCharge() * Acts::UnitConstants::e,
0033                  from.getMass() * Acts::UnitConstants::GeV);
0034
0035   // TODO do we have that in EDM4hep?
0036   // particle.setProcess(static_cast<ActsFatras::ProcessType>(data.process));
0037
0038   to.initialState().setPosition4(from.getVertex()[0] * Acts::UnitConstants::mm,
0039                                  from.getVertex()[1] * Acts::UnitConstants::mm,
0040                                  from.getVertex()[2] * Acts::UnitConstants::mm,
0041                                  from.getTime() * Acts::UnitConstants::ns);
0042
0043   // Only used for direction; normalization/units do not matter
0044   Acts::Vector3 momentum = {from.getMomentum()[0], from.getMomentum()[1],
0045                             from.getMomentum()[2]};
0046   to.initialState().setDirection(momentum.normalized());
0047
0048   to.initialState().setAbsoluteMomentum(momentum.norm() * 1_GeV);
0049
0050   return to;
0051 }
0052
0053 void EDM4hepUtil::writeParticle(const SimParticle& from,
0054                                 edm4hep::MutableMCParticle to) {
0055   // TODO what about particleId?
0056
0057   to.setPDG(from.pdg());
0058   to.setCharge(from.charge() / Acts::UnitConstants::e);
0059   to.setMass(from.mass() / Acts::UnitConstants::GeV);
0060   to.setVertex({from.position().x(), from.position().y(), from.position().z()});
0061   to.setMomentum({static_cast<float>(from.fourMomentum().x()),
0062                   static_cast<float>(from.fourMomentum().y()),
0063                   static_cast<float>(from.fourMomentum().z())});
0064   to.setMomentumAtEndpoint(
0065       {static_cast<float>(from.finalState().fourMomentum().x()),
0066        static_cast<float>(from.finalState().fourMomentum().y()),
0067        static_cast<float>(from.finalState().fourMomentum().z())});
0068 }
0069
0070 ActsFatras::Hit EDM4hepUtil::readSimHit(
0071     const edm4hep::SimTrackerHit& from, const MapParticleIdFrom& particleMapper,
0072     const MapGeometryIdFrom& geometryMapper) {
0073   auto particle = ActsPlugins::EDM4hepUtil::getParticle(from);
0074   ActsFatras::Barcode particleId = particleMapper(particle);
0075
0076   const auto mass = particle.getMass() * 1_GeV;
0077   const Acts::Vector3 momentum{
0078       from.getMomentum().x * 1_GeV,
0079       from.getMomentum().y * 1_GeV,
0080       from.getMomentum().z * 1_GeV,
0081   };
0082   const auto energy = std::hypot(momentum.norm(), mass);
0083
0084   Acts::Vector4 pos4{
0085       from.getPosition().x * 1_mm,
0086       from.getPosition().y * 1_mm,
0087       from.getPosition().z * 1_mm,
0088       from.getTime() * 1_ns,
0089   };
0090
0091   Acts::Vector4 mom4{
0092       momentum.x(),
0093       momentum.y(),
0094       momentum.z(),
0095       energy,
0096   };
0097
0098   Acts::Vector4 delta4 = Acts::Vector4::Zero();
0099   delta4[Acts::eEnergy] = -from.getEDep() * Acts::UnitConstants::GeV;
0100
0101   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0102
0103   // Can extract from time, but we need a complete picture of the trajectory
0104   // first
0105   std::int32_t index = -1;
0106
0107   return ActsFatras::Hit(geometryId, particleId, pos4, mom4, mom4 + delta4,
0108                          index);
0109 }
0110
0111 void EDM4hepUtil::writeSimHit(const ActsFatras::Hit& from,
0112                               edm4hep::MutableSimTrackerHit to,
0113                               const MapParticleIdTo& particleMapper,
0114                               const MapGeometryIdTo& geometryMapper) {
0115   const Acts::Vector4& globalPos4 = from.fourPosition();
0116   const Acts::Vector4& momentum4Before = from.momentum4Before();
0117   const auto delta4 = from.momentum4After() - momentum4Before;
0118
0119   if (particleMapper) {
0120     ActsPlugins::EDM4hepUtil::setParticle(to,
0121                                           particleMapper(from.particleId()));
0122   }
0123
0124   if (geometryMapper) {
0125     // TODO what about the digitization?
0126     to.setCellID(geometryMapper(from.geometryId()));
0127   }
0128
0129   to.setTime(globalPos4[Acts::eTime] / Acts::UnitConstants::ns);
0130
0131   to.setPosition({
0132       globalPos4[Acts::ePos0] / Acts::UnitConstants::mm,
0133       globalPos4[Acts::ePos1] / Acts::UnitConstants::mm,
0134       globalPos4[Acts::ePos2] / Acts::UnitConstants::mm,
0135   });
0136
0137   to.setMomentum({
0138       static_cast<float>(momentum4Before[Acts::eMom0] /
0139                          Acts::UnitConstants::GeV),
0140       static_cast<float>(momentum4Before[Acts::eMom1] /
0141                          Acts::UnitConstants::GeV),
0142       static_cast<float>(momentum4Before[Acts::eMom2] /
0143                          Acts::UnitConstants::GeV),
0144   });
0145
0146   to.setEDep(-delta4[Acts::eEnergy] / Acts::UnitConstants::GeV);
0147 }
0148
0149 VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement(
0150     MeasurementContainer& container, const edm4hep::TrackerHitPlane& from,
0151     const edm4hep::TrackerHit3DCollection* /*fromClusters*/,
0152     Cluster* /*toCluster*/, const MapGeometryIdFrom& geometryMapper) {
0153   // no need for digitization as we only want to identify the sensor
0154   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0155
0156   auto pos = from.getPosition();
0157   auto cov = from.getCovMatrix();
0158
0159   DigitizedParameters dParameters;
0160
0161   dParameters.indices.push_back(Acts::eBoundLoc0);
0162   dParameters.values.push_back(pos.x);
0163   dParameters.variances.push_back(cov[0]);
0164
0165   // TODO cut this out for 1D
0166   dParameters.indices.push_back(Acts::eBoundLoc1);
0167   dParameters.values.push_back(pos.y);
0168   dParameters.variances.push_back(cov[2]);
0169
0170   dParameters.indices.push_back(Acts::eBoundTime);
0171   dParameters.values.push_back(pos.z);
0172   dParameters.variances.push_back(cov[5]);
0173
0174   auto to = createMeasurement(container, geometryId, dParameters);
0175
0176   // @TODO: Figure out if cell information is accessible
0177
0178   return to;
0179 }
0180
0181 void EDM4hepUtil::writeMeasurement(
0182     const ConstVariableBoundMeasurementProxy& from,
0183     edm4hep::MutableTrackerHitPlane to, const Cluster* /*fromCluster*/,
0184     edm4hep::TrackerHit3DCollection& /*toClusters*/,
0185     const MapGeometryIdTo& geometryMapper) {
0186   Acts::GeometryIdentifier geoId = from.geometryId();
0187
0188   if (geometryMapper) {
0189     // no need for digitization as we only want to identify the sensor
0190     to.setCellID(geometryMapper(geoId));
0191   }
0192
0193   const auto& parameters = from.fullParameters();
0194   const auto& covariance = from.fullCovariance();
0195
0196   to.setTime(parameters[Acts::eBoundTime] / Acts::UnitConstants::ns);
0197
0198   to.setType(ActsPlugins::EDM4hepUtil::EDM4HEP_ACTS_POSITION_TYPE);
0199   // TODO set uv (which are in global spherical coordinates with r=1)
0200   to.setPosition({parameters[Acts::eBoundLoc0], parameters[Acts::eBoundLoc1],
0201                   parameters[Acts::eBoundTime]});
0202
0203   to.setCovMatrix({
0204       static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0205       static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc0)),
0206       static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0207       0,
0208       0,
0209       0,
0210   });
0211
0212   // @TODO: Check if we can write cell info
0213 }
0214
0215 void EDM4hepUtil::writeTrajectory(
0216     const Acts::GeometryContext& gctx, double Bz, const Trajectories& from,
0217     edm4hep::MutableTrack to, std::size_t fromIndex,
0218     const Acts::ParticleHypothesis& particleHypothesis,
0219     const IndexMultimap<ActsFatras::Barcode>& hitParticlesMap) {
0220   const auto& multiTrajectory = from.multiTrajectory();
0221   auto trajectoryState =
0222       Acts::MultiTrajectoryHelpers::trajectoryState(multiTrajectory, fromIndex);
0223
0224   std::vector<ParticleHitCount> particleHitCount;
0225   identifyContributingParticles(hitParticlesMap, from, fromIndex,
0226                                 particleHitCount);
0227   // TODO use particles
0228
0229   // TODO write track params
0230   // auto trackParameters = from.trackParameters(fromIndex);
0231
0232   to.setChi2(trajectoryState.chi2Sum / trajectoryState.NDF);
0233   to.setNdf(trajectoryState.NDF);
0234
0235   multiTrajectory.visitBackwards(fromIndex, [&](const auto& state) {
0236     // we only fill the track states with non-outlier measurement
0237     auto typeFlags = state.typeFlags();
0238     if (!typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
0239       return true;
0240     }
0241
0242     edm4hep::TrackState trackState;
0243
0244     Acts::BoundTrackParameters parObj{state.referenceSurface().getSharedPtr(),
0245                                       state.parameters(), state.covariance(),
0246                                       particleHypothesis};
0247
0248     // Convert to LCIO track parametrization expected by EDM4hep
0249     // This will create an ad-hoc perigee surface if the input parameters are
0250     // not bound on a perigee surface already
0251     ActsPlugins::EDM4hepUtil::detail::Parameters converted =
0252         ActsPlugins::EDM4hepUtil::detail::convertTrackParametersToEdm4hep(
0253             gctx, Bz, parObj);
0254
0255     trackState.D0 = converted.values[0];
0256     trackState.Z0 = converted.values[1];
0257     trackState.phi = converted.values[2];
0258     trackState.tanLambda = converted.values[3];
0259     trackState.omega = converted.values[4];
0260     trackState.time = converted.values[5];
0261
0262     // Converted parameters are relative to an ad-hoc perigee surface created at
0263     // the hit location
0264     auto center = converted.surface->center(gctx);
0265     trackState.referencePoint.x = center.x();
0266     trackState.referencePoint.y = center.y();
0267     trackState.referencePoint.z = center.z();
0268
0269     if (converted.covariance) {
0270       auto c = [&](std::size_t row, std::size_t col) {
0271         return static_cast<float>(converted.covariance.value()(row, col));
0272       };
0273
0274       // clang-format off
0275       trackState.covMatrix = {c(0, 0),
0276                               c(1, 0), c(1, 1),
0277                               c(2, 0), c(2, 1), c(2, 2),
0278                               c(3, 0), c(3, 1), c(3, 2), c(3, 3),
0279                               c(4, 0), c(4, 1), c(4, 2), c(4, 3), c(4, 4),
0280                               c(5, 0), c(5, 1), c(5, 2), c(5, 3), c(5, 4), c(5, 5)};
0281       // clang-format on
0282     }
0283
0284     to.addToTrackStates(trackState);
0285
0286     return true;
0287   });
0288 }
0289
0290 }  // namespace ActsExamples