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