EDM4hep/src/EDM4hepUtil.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2022 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 http://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/Charge.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
0018 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
0019 #include "ActsExamples/EventData/Index.hpp"
0020 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0021 #include "ActsExamples/EventData/SimHit.hpp"
0022 #include "ActsExamples/Validation/TrackClassification.hpp"
0023
0024 #include "edm4hep/TrackState.h"
0025
0026 using namespace Acts::UnitLiterals;
0027
0028 namespace ActsExamples {
0029
0030 ActsFatras::Particle EDM4hepUtil::readParticle(
0031     const edm4hep::MCParticle& from, const MapParticleIdFrom& particleMapper) {
0032   ActsFatras::Barcode particleId = particleMapper(from);
0033
0034   ActsFatras::Particle to(particleId,
0035                           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.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.setDirection(momentum.normalized());
0051
0052   to.setAbsoluteMomentum(momentum.norm() * 1_GeV);
0053
0054   return to;
0055 }
0056
0057 void EDM4hepUtil::writeParticle(const ActsFatras::Particle& 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 }
0069
0070 ActsFatras::Hit EDM4hepUtil::readSimHit(
0071     const edm4hep::SimTrackerHit& from, const MapParticleIdFrom& particleMapper,
0072     const MapGeometryIdFrom& geometryMapper) {
0073   ActsFatras::Barcode particleId = particleMapper(from.getMCParticle());
0074
0075   const auto mass = from.getMCParticle().getMass() * 1_GeV;
0076   const Acts::Vector3 momentum{
0077       from.getMomentum().x * 1_GeV,
0078       from.getMomentum().y * 1_GeV,
0079       from.getMomentum().z * 1_GeV,
0080   };
0081   const auto energy = std::hypot(momentum.norm(), mass);
0082
0083   Acts::Vector4 pos4{
0084       from.getPosition().x * 1_mm,
0085       from.getPosition().y * 1_mm,
0086       from.getPosition().z * 1_mm,
0087       from.getTime() * 1_ns,
0088   };
0089
0090   Acts::Vector4 mom4{
0091       momentum.x(),
0092       momentum.y(),
0093       momentum.z(),
0094       energy,
0095   };
0096
0097   Acts::Vector4 delta4 = Acts::Vector4::Zero();
0098   delta4[Acts::eEnergy] = -from.getEDep() * Acts::UnitConstants::GeV;
0099
0100   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0101
0102   // Can extract from time, but we need a complete picture of the trajectory
0103   // first
0104   int32_t index = -1;
0105
0106   return ActsFatras::Hit(geometryId, particleId, pos4, mom4, mom4 + delta4,
0107                          index);
0108 }
0109
0110 void EDM4hepUtil::writeSimHit(const ActsFatras::Hit& from,
0111                               edm4hep::MutableSimTrackerHit to,
0112                               const MapParticleIdTo& particleMapper,
0113                               const MapGeometryIdTo& geometryMapper) {
0114   const Acts::Vector4& globalPos4 = from.fourPosition();
0115   const Acts::Vector4& momentum4Before = from.momentum4Before();
0116   const auto delta4 = from.momentum4After() - momentum4Before;
0117
0118   if (particleMapper) {
0119     to.setMCParticle(particleMapper(from.particleId()));
0120   }
0121
0122   if (geometryMapper) {
0123     // TODO what about the digitization?
0124     to.setCellID(geometryMapper(from.geometryId()));
0125   }
0126
0127   to.setTime(globalPos4[Acts::eTime] / Acts::UnitConstants::ns);
0128
0129   to.setPosition({
0130       globalPos4[Acts::ePos0] / Acts::UnitConstants::mm,
0131       globalPos4[Acts::ePos1] / Acts::UnitConstants::mm,
0132       globalPos4[Acts::ePos2] / Acts::UnitConstants::mm,
0133   });
0134
0135   to.setMomentum({
0136       static_cast<float>(momentum4Before[Acts::eMom0] /
0137                          Acts::UnitConstants::GeV),
0138       static_cast<float>(momentum4Before[Acts::eMom1] /
0139                          Acts::UnitConstants::GeV),
0140       static_cast<float>(momentum4Before[Acts::eMom2] /
0141                          Acts::UnitConstants::GeV),
0142   });
0143
0144   to.setEDep(-delta4[Acts::eEnergy] / Acts::UnitConstants::GeV);
0145 }
0146
0147 Measurement EDM4hepUtil::readMeasurement(
0148     const edm4hep::TrackerHitPlane& from,
0149     const edm4hep::TrackerHitCollection* fromClusters, Cluster* toCluster,
0150     const MapGeometryIdFrom& geometryMapper) {
0151   // no need for digitization as we only want to identify the sensor
0152   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
0153
0154   IndexSourceLink sourceLink{
0155       geometryId, static_cast<Index>(podioObjectIDToInteger(from.id()))};
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(dParameters, sourceLink);
0176
0177   if (fromClusters != nullptr) {
0178     for (const auto objectId : from.getRawHits()) {
0179       const auto& c = fromClusters->at(objectId.index);
0180
0181       // TODO get EDM4hep fixed
0182       // misusing some fields to store ACTS specific information
0183       // don't ask ...
0184       ActsFatras::Segmentizer::Bin2D bin{
0185           static_cast<unsigned int>(c.getType()),
0186           static_cast<unsigned int>(c.getQuality())};
0187       ActsFatras::Segmentizer::Segment2D path2D{
0188           {Acts::Vector2::Zero(), Acts::Vector2::Zero()}};
0189       double activation = c.getTime();
0190       ActsFatras::Segmentizer::ChannelSegment cell{bin, path2D, activation};
0191
0192       toCluster->channels.push_back(cell);
0193     }
0194   }
0195
0196   return to;
0197 }
0198
0199 void EDM4hepUtil::writeMeasurement(const Measurement& from,
0200                                    edm4hep::MutableTrackerHitPlane to,
0201                                    const Cluster* fromCluster,
0202                                    edm4hep::TrackerHitCollection& toClusters,
0203                                    const MapGeometryIdTo& geometryMapper) {
0204   std::visit(
0205       [&](const auto& m) {
0206         Acts::GeometryIdentifier geoId =
0207             m.sourceLink().template get<IndexSourceLink>().geometryId();
0208
0209         if (geometryMapper) {
0210           // no need for digitization as we only want to identify the sensor
0211           to.setCellID(geometryMapper(geoId));
0212         }
0213
0214         auto parameters = (m.expander() * m.parameters()).eval();
0215
0216         to.setTime(parameters[Acts::eBoundTime] / Acts::UnitConstants::ns);
0217
0218         to.setType(Acts::EDM4hepUtil::EDM4HEP_ACTS_POSITION_TYPE);
0219         // TODO set uv (which are in global spherical coordinates with r=1)
0220         to.setPosition({parameters[Acts::eBoundLoc0],
0221                         parameters[Acts::eBoundLoc1],
0222                         parameters[Acts::eBoundTime]});
0223
0224         auto covariance =
0225             (m.expander() * m.covariance() * m.expander().transpose()).eval();
0226         to.setCovMatrix({
0227             static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0228             static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc0)),
0229             static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0230             0,
0231             0,
0232             0,
0233         });
0234
0235         if (fromCluster) {
0236           for (const auto& c : fromCluster->channels) {
0237             auto toChannel = toClusters.create();
0238             to.addToRawHits(toChannel.getObjectID());
0239
0240             // TODO digitization channel
0241
0242             // TODO get EDM4hep fixed
0243             // misusing some fields to store ACTS specific information
0244             // don't ask ...
0245             toChannel.setType(c.bin[0]);
0246             toChannel.setQuality(c.bin[1]);
0247             toChannel.setTime(c.activation);
0248           }
0249         }
0250       },
0251       from);
0252 }
0253
0254 void EDM4hepUtil::writeTrajectory(
0255     const Acts::GeometryContext& gctx, double Bz, const Trajectories& from,
0256     edm4hep::MutableTrack to, std::size_t fromIndex,
0257     const Acts::ParticleHypothesis& particleHypothesis,
0258     const IndexMultimap<ActsFatras::Barcode>& hitParticlesMap) {
0259   const auto& multiTrajectory = from.multiTrajectory();
0260   auto trajectoryState =
0261       Acts::MultiTrajectoryHelpers::trajectoryState(multiTrajectory, fromIndex);
0262
0263   std::vector<ParticleHitCount> particleHitCount;
0264   identifyContributingParticles(hitParticlesMap, from, fromIndex,
0265                                 particleHitCount);
0266   // TODO use particles
0267
0268   // TODO write track params
0269   // auto trackParameters = from.trackParameters(fromIndex);
0270
0271   to.setChi2(trajectoryState.chi2Sum / trajectoryState.NDF);
0272   to.setNdf(trajectoryState.NDF);
0273
0274   multiTrajectory.visitBackwards(fromIndex, [&](const auto& state) {
0275     // we only fill the track states with non-outlier measurement
0276     auto typeFlags = state.typeFlags();
0277     if (!typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
0278       return true;
0279     }
0280
0281     edm4hep::TrackState trackState;
0282
0283     Acts::BoundTrackParameters parObj{state.referenceSurface().getSharedPtr(),
0284                                       state.parameters(), state.covariance(),
0285                                       particleHypothesis};
0286
0287     // Convert to LCIO track parametrization expected by EDM4hep
0288     // This will create an ad-hoc perigee surface if the input parameters are
0289     // not bound on a perigee surface already
0290     Acts::EDM4hepUtil::detail::Parameters converted =
0291         Acts::EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz,
0292                                                                    parObj);
0293
0294     trackState.D0 = converted.values[0];
0295     trackState.Z0 = converted.values[1];
0296     trackState.phi = converted.values[2];
0297     trackState.tanLambda = converted.values[3];
0298     trackState.omega = converted.values[4];
0299     trackState.time = converted.values[5];
0300
0301     // Converted parameters are relative to an ad-hoc perigee surface created at
0302     // the hit location
0303     auto center = converted.surface->center(gctx);
0304     trackState.referencePoint.x = center.x();
0305     trackState.referencePoint.y = center.y();
0306     trackState.referencePoint.z = center.z();
0307
0308     if (converted.covariance) {
0309       const auto& c = converted.covariance.value();
0310
0311       trackState.covMatrix = {
0312           static_cast<float>(c(0, 0)), static_cast<float>(c(1, 0)),
0313           static_cast<float>(c(1, 1)), static_cast<float>(c(2, 0)),
0314           static_cast<float>(c(2, 1)), static_cast<float>(c(2, 2)),
0315           static_cast<float>(c(3, 0)), static_cast<float>(c(3, 1)),
0316           static_cast<float>(c(3, 2)), static_cast<float>(c(3, 3)),
0317           static_cast<float>(c(4, 0)), static_cast<float>(c(4, 1)),
0318           static_cast<float>(c(4, 2)), static_cast<float>(c(4, 3)),
0319           static_cast<float>(c(4, 4))};
0320     }
0321
0322     to.addToTrackStates(trackState);
0323
0324     return true;
0325   });
0326 }
0327
0328 }  // namespace ActsExamples