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