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& initialState() const { return m_initial; }
0060   const SimParticleState& finalState() const { return m_final; }
0061
0062   SimParticleState& initialState() { return m_initial; }
0063   SimParticleState& finalState() { 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(initialState().withParticleId(particleId),
0072                        finalState().withParticleId(particleId));
0073   }
0074
0075   /// Set the process type that generated this particle.
0076   SimParticle& setProcess(ActsFatras::ProcessType proc) {
0077     initialState().setProcess(proc);
0078     finalState().setProcess(proc);
0079     return *this;
0080   }
0081   /// Set the pdg.
0082   SimParticle& setPdg(Acts::PdgParticle pdg) {
0083     initialState().setPdg(pdg);
0084     finalState().setPdg(pdg);
0085     return *this;
0086   }
0087   /// Set the charge.
0088   SimParticle& setCharge(double charge) {
0089     initialState().setCharge(charge);
0090     finalState().setCharge(charge);
0091     return *this;
0092   }
0093   /// Set the mass.
0094   SimParticle& setMass(double mass) {
0095     initialState().setMass(mass);
0096     finalState().setMass(mass);
0097     return *this;
0098   }
0099   /// Set the particle ID.
0100   SimParticle& setParticleId(SimBarcode barcode) {
0101     initialState().setParticleId(barcode);
0102     finalState().setParticleId(barcode);
0103     return *this;
0104   }
0105
0106   /// Particle identifier within an event.
0107   SimBarcode particleId() const { return initialState().particleId(); }
0108   /// Which type of process generated this particle.
0109   ActsFatras::ProcessType process() const { return initialState().process(); }
0110   /// PDG particle number that identifies the type.
0111   Acts::PdgParticle pdg() const { return initialState().pdg(); }
0112   /// Absolute PDG particle number that identifies the type.
0113   Acts::PdgParticle absolutePdg() const { return initialState().absolutePdg(); }
0114   /// Particle charge.
0115   double charge() const { return initialState().charge(); }
0116   /// Particle absolute charge.
0117   double absoluteCharge() const { return initialState().absoluteCharge(); }
0118   /// Particle mass.
0119   double mass() const { return initialState().mass(); }
0120
0121   /// Check if this is a secondary particle.
0122   bool isSecondary() const { return initialState().isSecondary(); }
0123
0124   /// Particle hypothesis.
0125   Acts::ParticleHypothesis hypothesis() const {
0126     return initialState().hypothesis();
0127   }
0128   /// Particl qOverP.
0129   double qOverP() const { return initialState().qOverP(); }
0130
0131   /// Space-time position four-vector.
0132   const Acts::Vector4& fourPosition() const {
0133     return initialState().fourPosition();
0134   }
0135   /// Three-position, i.e. spatial coordinates without the time.
0136   auto position() const { return initialState().position(); }
0137   /// Time coordinate.
0138   double time() const { return initialState().time(); }
0139   /// Energy-momentum four-vector.
0140   Acts::Vector4 fourMomentum() const { return initialState().fourMomentum(); }
0141   /// Unit three-direction, i.e. the normalized momentum three-vector.
0142   const Acts::Vector3& direction() const { return initialState().direction(); }
0143   /// Polar angle.
0144   double theta() const { return initialState().theta(); }
0145   /// Azimuthal angle.
0146   double phi() const { return initialState().phi(); }
0147   /// Absolute momentum in the x-y plane.
0148   double transverseMomentum() const {
0149     return initialState().transverseMomentum();
0150   }
0151   /// Absolute momentum.
0152   double absoluteMomentum() const { return initialState().absoluteMomentum(); }
0153   /// Absolute momentum.
0154   Acts::Vector3 momentum() const { return initialState().momentum(); }
0155   /// Total energy, i.e. norm of the four-momentum.
0156   double energy() const { return initialState().energy(); }
0157
0158   /// Energy loss over the particles lifetime or simulation time.
0159   double energyLoss() const {
0160     return initialState().energy() - finalState().energy();
0161   }
0162
0163   /// Accumulated path within material measured in radiation lengths.
0164   double pathInX0() const { return finalState().pathInX0(); }
0165   /// Accumulated path within material measured in interaction lengths.
0166   double pathInL0() const { return finalState().pathInL0(); }
0167
0168   /// Number of hits.
0169   std::uint32_t numberOfHits() const { return finalState().numberOfHits(); }
0170
0171   /// Particle outcome.
0172   ActsFatras::ParticleOutcome outcome() const { return finalState().outcome(); }
0173
0174  private:
0175   SimParticleState m_initial;
0176   SimParticleState m_final;
0177 };
0178
0179 std::ostream& operator<<(std::ostream& os, const SimParticle& particle);
0180
0181 namespace detail {
0182 struct CompareParticleId {
0183   using is_transparent = void;
0184   bool operator()(const SimParticleState& lhs,
0185                   const SimParticleState& rhs) const {
0186     return lhs.particleId() < rhs.particleId();
0187   }
0188   bool operator()(const SimParticle& lhs, const SimParticle& rhs) const {
0189     return lhs.particleId() < rhs.particleId();
0190   }
0191   bool operator()(SimBarcode lhs, const SimParticleState& rhs) const {
0192     return lhs < rhs.particleId();
0193   }
0194   bool operator()(SimBarcode lhs, const SimParticle& rhs) const {
0195     return lhs < rhs.particleId();
0196   }
0197   bool operator()(const SimParticleState& lhs, SimBarcode rhs) const {
0198     return lhs.particleId() < rhs;
0199   }
0200   bool operator()(const SimParticle& lhs, SimBarcode rhs) const {
0201     return lhs.particleId() < rhs;
0202   }
0203 };
0204 struct PrimaryVertexIdGetter {
0205   SimBarcode operator()(const SimParticleState& particle) const {
0206     return SimBarcode().withVertexPrimary(
0207         particle.particleId().vertexPrimary());
0208   }
0209   SimBarcode operator()(const SimParticle& particle) const {
0210     return SimBarcode().withVertexPrimary(
0211         particle.particleId().vertexPrimary());
0212   }
0213 };
0214 struct SecondaryVertexIdGetter {
0215   SimBarcode operator()(const SimParticleState& particle) const {
0216     return SimBarcode()
0217         .withVertexPrimary(particle.particleId().vertexPrimary())
0218         .withVertexSecondary(particle.particleId().vertexSecondary());
0219   }
0220   SimBarcode operator()(const SimParticle& particle) const {
0221     return SimBarcode()
0222         .withVertexPrimary(particle.particleId().vertexPrimary())
0223         .withVertexSecondary(particle.particleId().vertexSecondary());
0224   }
0225 };
0226 struct VertexIdGetter {
0227   SimBarcode operator()(const SimParticleState& particle) const {
0228     return particle.particleId().vertexId();
0229   }
0230   SimBarcode operator()(const SimParticle& particle) const {
0231     return particle.particleId().vertexId();
0232   }
0233 };
0234 }  // namespace detail
0235
0236 using SimParticleStateContainer =
0237     ::boost::container::flat_set<SimParticleState, detail::CompareParticleId>;
0238
0239 /// Store particles ordered by particle identifier.
0240 using SimParticleContainer =
0241     ::boost::container::flat_set<SimParticle, detail::CompareParticleId>;
0242
0243 /// Iterate over groups of particles belonging to the same primary vertex.
0244 inline GroupBy<SimParticleContainer::const_iterator,
0245                detail::PrimaryVertexIdGetter>
0246 groupByPrimaryVertex(const SimParticleContainer& container) {
0247   return makeGroupBy(container, detail::PrimaryVertexIdGetter());
0248 }
0249
0250 /// Iterate over groups of particles belonging to the same secondary vertex.
0251 ///
0252 /// For each primary vertex, this yields one group of particles belonging
0253 /// directly to the secondary vertex and a group for each secondary vertex.
0254 inline GroupBy<SimParticleContainer::const_iterator,
0255                detail::SecondaryVertexIdGetter>
0256 groupBySecondaryVertex(const SimParticleContainer& container) {
0257   return makeGroupBy(container, detail::SecondaryVertexIdGetter());
0258 }
0259
0260 /// Iterate over groups of particles belonging to the same vertex.
0261 ///
0262 /// For each vertex, this yields one group of particles belonging
0263 /// directly to the vertex and a group for each secondary vertex.
0264 inline GroupBy<SimParticleContainer::const_iterator, detail::VertexIdGetter>
0265 groupByVertexId(const SimParticleContainer& container) {
0266   return makeGroupBy(container, detail::VertexIdGetter());
0267 }
0268
0269 }  // namespace ActsExamples