ActsExamples/EventData/SimParticle.hpp

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 #pragma once
0010
0011 #include "ActsExamples/Utilities/GroupBy.hpp"
0012 #include "ActsFatras/EventData/Particle.hpp"
0013 #include "ActsFatras/EventData/ParticleOutcome.hpp"
0014
0015 #include <boost/container/flat_set.hpp>
0016
0017 namespace ActsExamples {
0018
0019 using SimBarcode = ::ActsFatras::Barcode;
0020 using SimBarcodeContainer = ::boost::container::flat_set<SimBarcode>;
0021
0022 using SimParticleState = ::ActsFatras::Particle;
0023
0024 class SimParticle final {
0025  public:
0026   /// Construct a default particle with invalid identity.
0027   SimParticle() = default;
0028
0029   /// Construct a particle at rest with explicit mass and charge.
0030   ///
0031   /// @param particleId Particle identifier within an event
0032   /// @param pdg PDG id
0033   /// @param charge Particle charge in native units
0034   /// @param mass Particle mass in native units
0035   ///
0036   /// @warning It is the users responsibility that charge and mass match
0037   ///          the PDG particle number.
0038   SimParticle(SimBarcode particleId, Acts::PdgParticle pdg, double charge,
0039               double mass)
0040       : m_initial(particleId, pdg, charge, mass),
0041         m_final(particleId, pdg, charge, mass) {}
0042
0043   /// Construct a particle at rest from a PDG particle number.
0044   ///
0045   /// @param particleId Particle identifier within an event
0046   /// @param pdg PDG particle number
0047   ///
0048   /// Charge and mass are retrieved from the particle data table.
0049   SimParticle(SimBarcode particleId, Acts::PdgParticle pdg)
0050       : m_initial(particleId, pdg), m_final(particleId, pdg) {}
0051
0052   SimParticle(const SimParticleState& initial, const SimParticleState& final)
0053       : m_initial(initial), m_final(final) {
0054     if (m_initial.particleId() != m_final.particleId()) {
0055       throw std::invalid_argument("Particle id mismatch");
0056     }
0057   }
0058
0059   const SimParticleState& initial() const { return m_initial; }
0060   const SimParticleState& final() const { return m_final; }
0061
0062   SimParticleState& initial() { return m_initial; }
0063   SimParticleState& final() { return m_final; }
0064
0065   /// Construct a new particle with a new identifier but same kinematics.
0066   ///
0067   /// @note This is intentionally not a regular setter. The particle id
0068   ///       is used to identify the whole particle. Setting it on an existing
0069   ///       particle is usually a mistake.
0070   SimParticle withParticleId(SimBarcode particleId) const {
0071     return SimParticle(initial().withParticleId(particleId),
0072                        final().withParticleId(particleId));
0073   }
0074
0075   /// Set the process type that generated this particle.
0076   SimParticle& setProcess(ActsFatras::ProcessType proc) {
0077     initial().setProcess(proc);
0078     final().setProcess(proc);
0079     return *this;
0080   }
0081   /// Set the pdg.
0082   SimParticle& setPdg(Acts::PdgParticle pdg) {
0083     initial().setPdg(pdg);
0084     final().setPdg(pdg);
0085     return *this;
0086   }
0087   /// Set the charge.
0088   SimParticle& setCharge(double charge) {
0089     initial().setCharge(charge);
0090     final().setCharge(charge);
0091     return *this;
0092   }
0093   /// Set the mass.
0094   SimParticle& setMass(double mass) {
0095     initial().setMass(mass);
0096     final().setMass(mass);
0097     return *this;
0098   }
0099   /// Set the particle ID.
0100   SimParticle& setParticleId(SimBarcode barcode) {
0101     initial().setParticleId(barcode);
0102     final().setParticleId(barcode);
0103     return *this;
0104   }
0105
0106   /// Particle identifier within an event.
0107   SimBarcode particleId() const { return initial().particleId(); }
0108   /// Which type of process generated this particle.
0109   ActsFatras::ProcessType process() const { return initial().process(); }
0110   /// PDG particle number that identifies the type.
0111   Acts::PdgParticle pdg() const { return initial().pdg(); }
0112   /// Absolute PDG particle number that identifies the type.
0113   Acts::PdgParticle absolutePdg() const { return initial().absolutePdg(); }
0114   /// Particle charge.
0115   double charge() const { return initial().charge(); }
0116   /// Particle absolute charge.
0117   double absoluteCharge() const { return initial().absoluteCharge(); }
0118   /// Particle mass.
0119   double mass() const { return initial().mass(); }
0120
0121   /// Check if this is a secondary particle.
0122   bool isSecondary() const { return initial().isSecondary(); }
0123
0124   /// Particle hypothesis.
0125   Acts::ParticleHypothesis hypothesis() const { return initial().hypothesis(); }
0126   /// Particl qOverP.
0127   double qOverP() const { return initial().qOverP(); }
0128
0129   /// Space-time position four-vector.
0130   const Acts::Vector4& fourPosition() const { return initial().fourPosition(); }
0131   /// Three-position, i.e. spatial coordinates without the time.
0132   auto position() const { return initial().position(); }
0133   /// Time coordinate.
0134   double time() const { return initial().time(); }
0135   /// Energy-momentum four-vector.
0136   Acts::Vector4 fourMomentum() const { return initial().fourMomentum(); }
0137   /// Unit three-direction, i.e. the normalized momentum three-vector.
0138   const Acts::Vector3& direction() const { return initial().direction(); }
0139   /// Polar angle.
0140   double theta() const { return initial().theta(); }
0141   /// Azimuthal angle.
0142   double phi() const { return initial().phi(); }
0143   /// Absolute momentum in the x-y plane.
0144   double transverseMomentum() const { return initial().transverseMomentum(); }
0145   /// Absolute momentum.
0146   double absoluteMomentum() const { return initial().absoluteMomentum(); }
0147   /// Absolute momentum.
0148   Acts::Vector3 momentum() const { return initial().momentum(); }
0149   /// Total energy, i.e. norm of the four-momentum.
0150   double energy() const { return initial().energy(); }
0151
0152   /// Energy loss over the particles lifetime or simulation time.
0153   double energyLoss() const { return initial().energy() - final().energy(); }
0154
0155   /// Accumulated path within material measured in radiation lengths.
0156   double pathInX0() const { return final().pathInX0(); }
0157   /// Accumulated path within material measured in interaction lengths.
0158   double pathInL0() const { return final().pathInL0(); }
0159
0160   /// Number of hits.
0161   std::uint32_t numberOfHits() const { return final().numberOfHits(); }
0162
0163   /// Particle outcome.
0164   ActsFatras::ParticleOutcome outcome() const { return final().outcome(); }
0165
0166  private:
0167   SimParticleState m_initial;
0168   SimParticleState m_final;
0169 };
0170
0171 std::ostream& operator<<(std::ostream& os, const SimParticle& particle);
0172
0173 namespace detail {
0174 struct CompareParticleId {
0175   using is_transparent = void;
0176   bool operator()(const SimParticleState& lhs,
0177                   const SimParticleState& rhs) const {
0178     return lhs.particleId() < rhs.particleId();
0179   }
0180   bool operator()(const SimParticle& lhs, const SimParticle& rhs) const {
0181     return lhs.particleId() < rhs.particleId();
0182   }
0183   bool operator()(SimBarcode lhs, const SimParticleState& rhs) const {
0184     return lhs < rhs.particleId();
0185   }
0186   bool operator()(SimBarcode lhs, const SimParticle& rhs) const {
0187     return lhs < rhs.particleId();
0188   }
0189   bool operator()(const SimParticleState& lhs, SimBarcode rhs) const {
0190     return lhs.particleId() < rhs;
0191   }
0192   bool operator()(const SimParticle& lhs, SimBarcode rhs) const {
0193     return lhs.particleId() < rhs;
0194   }
0195 };
0196 struct PrimaryVertexIdGetter {
0197   SimBarcode operator()(const SimParticleState& particle) const {
0198     return SimBarcode(0u).setVertexPrimary(
0199         particle.particleId().vertexPrimary());
0200   }
0201   SimBarcode operator()(const SimParticle& particle) const {
0202     return SimBarcode(0u).setVertexPrimary(
0203         particle.particleId().vertexPrimary());
0204   }
0205 };
0206 struct SecondaryVertexIdGetter {
0207   SimBarcode operator()(const SimParticleState& particle) const {
0208     return SimBarcode(0u)
0209         .setVertexPrimary(particle.particleId().vertexPrimary())
0210         .setVertexSecondary(particle.particleId().vertexSecondary());
0211   }
0212   SimBarcode operator()(const SimParticle& particle) const {
0213     return SimBarcode(0u)
0214         .setVertexPrimary(particle.particleId().vertexPrimary())
0215         .setVertexSecondary(particle.particleId().vertexSecondary());
0216   }
0217 };
0218 }  // namespace detail
0219
0220 using SimParticleStateContainer =
0221     ::boost::container::flat_set<SimParticleState, detail::CompareParticleId>;
0222
0223 /// Store particles ordered by particle identifier.
0224 using SimParticleContainer =
0225     ::boost::container::flat_set<SimParticle, detail::CompareParticleId>;
0226
0227 /// Iterate over groups of particles belonging to the same primary vertex.
0228 inline GroupBy<SimParticleContainer::const_iterator,
0229                detail::PrimaryVertexIdGetter>
0230 groupByPrimaryVertex(const SimParticleContainer& container) {
0231   return makeGroupBy(container, detail::PrimaryVertexIdGetter());
0232 }
0233
0234 /// Iterate over groups of particles belonging to the same secondary vertex.
0235 ///
0236 /// For each primary vertex, this yields one group of particles belonging
0237 /// directly to the primary vertex (secondary vertex id 0) and a group for
0238 /// each secondary vertex.
0239 inline GroupBy<SimParticleContainer::const_iterator,
0240                detail::SecondaryVertexIdGetter>
0241 groupBySecondaryVertex(const SimParticleContainer& container) {
0242   return makeGroupBy(container, detail::SecondaryVertexIdGetter());
0243 }
0244
0245 }  // namespace ActsExamples