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