ActsFatras/EventData/Particle.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2024 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 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/PdgParticle.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/EventData/ParticleHypothesis.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Utilities/VectorHelpers.hpp"
0020 #include "ActsFatras/EventData/Barcode.hpp"
0021 #include "ActsFatras/EventData/ParticleOutcome.hpp"
0022 #include "ActsFatras/EventData/ProcessType.hpp"
0023
0024 #include <cmath>
0025 #include <iosfwd>
0026 #include <optional>
0027
0028 namespace ActsFatras {
0029
0030 /// Particle identity information and kinematic state.
0031 ///
0032 /// Also stores some simulation-specific properties.
0033 class Particle {
0034  public:
0035   using Scalar = Acts::ActsScalar;
0036   using Vector3 = Acts::ActsVector<3>;
0037   using Vector4 = Acts::ActsVector<4>;
0038
0039   /// Construct a default particle with invalid identity.
0040   Particle() = default;
0041   /// Construct a particle at rest with explicit mass and charge.
0042   ///
0043   /// @param particleId Particle identifier within an event
0044   /// @param pdg PDG id
0045   /// @param charge Particle charge in native units
0046   /// @param mass Particle mass in native units
0047   ///
0048   /// @warning It is the users responsibility that charge and mass match
0049   ///          the PDG particle number.
0050   Particle(Barcode particleId, Acts::PdgParticle pdg, Scalar charge,
0051            Scalar mass)
0052       : m_particleId(particleId), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
0053   /// Construct a particle at rest from a PDG particle number.
0054   ///
0055   /// @param particleId Particle identifier within an event
0056   /// @param pdg PDG particle number
0057   ///
0058   /// Charge and mass are retrieved from the particle data table.
0059   Particle(Barcode particleId, Acts::PdgParticle pdg);
0060   Particle(const Particle &) = default;
0061   Particle(Particle &&) = default;
0062   Particle &operator=(const Particle &) = default;
0063   Particle &operator=(Particle &&) = default;
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   Particle withParticleId(Barcode particleId) const {
0071     Particle p = *this;
0072     p.m_particleId = particleId;
0073     return p;
0074   }
0075
0076   /// Set the process type that generated this particle.
0077   Particle &setProcess(ProcessType proc) {
0078     m_process = proc;
0079     return *this;
0080   }
0081   /// Set the pdg.
0082   Particle setPdg(Acts::PdgParticle pdg) {
0083     m_pdg = pdg;
0084     return *this;
0085   }
0086   /// Set the charge.
0087   Particle setCharge(Scalar charge) {
0088     m_charge = charge;
0089     return *this;
0090   }
0091   /// Set the mass.
0092   Particle setMass(Scalar mass) {
0093     m_mass = mass;
0094     return *this;
0095   }
0096   /// Set the particle ID.
0097   Particle &setParticleId(Barcode barcode) {
0098     m_particleId = barcode;
0099     return *this;
0100   }
0101   /// Set the space-time position four-vector.
0102   Particle &setPosition4(const Vector4 &pos4) {
0103     m_position4 = pos4;
0104     return *this;
0105   }
0106   /// Set the space-time position four-vector from three-position and time.
0107   Particle &setPosition4(const Vector3 &position, Scalar time) {
0108     m_position4.segment<3>(Acts::ePos0) = position;
0109     m_position4[Acts::eTime] = time;
0110     return *this;
0111   }
0112   /// Set the space-time position four-vector from scalar components.
0113   Particle &setPosition4(Scalar x, Scalar y, Scalar z, Scalar time) {
0114     m_position4[Acts::ePos0] = x;
0115     m_position4[Acts::ePos1] = y;
0116     m_position4[Acts::ePos2] = z;
0117     m_position4[Acts::eTime] = time;
0118     return *this;
0119   }
0120   /// Set the direction three-vector
0121   Particle &setDirection(const Vector3 &direction) {
0122     m_direction = direction;
0123     m_direction.normalize();
0124     return *this;
0125   }
0126   /// Set the direction three-vector from scalar components.
0127   Particle &setDirection(Scalar dx, Scalar dy, Scalar dz) {
0128     m_direction[Acts::ePos0] = dx;
0129     m_direction[Acts::ePos1] = dy;
0130     m_direction[Acts::ePos2] = dz;
0131     m_direction.normalize();
0132     return *this;
0133   }
0134   /// Set the absolute momentum.
0135   Particle &setAbsoluteMomentum(Scalar absMomentum) {
0136     m_absMomentum = absMomentum;
0137     return *this;
0138   }
0139
0140   /// Change the energy by the given amount.
0141   ///
0142   /// Energy loss corresponds to a negative change. If the updated energy
0143   /// would result in an unphysical value, the particle is put to rest, i.e.
0144   /// its absolute momentum is set to zero.
0145   Particle &correctEnergy(Scalar delta) {
0146     const auto newEnergy = std::hypot(m_mass, m_absMomentum) + delta;
0147     if (newEnergy <= m_mass) {
0148       m_absMomentum = Scalar{0};
0149     } else {
0150       m_absMomentum = std::sqrt(newEnergy * newEnergy - m_mass * m_mass);
0151     }
0152     return *this;
0153   }
0154
0155   /// Particle identifier within an event.
0156   constexpr Barcode particleId() const { return m_particleId; }
0157   /// Which type of process generated this particle.
0158   constexpr ProcessType process() const { return m_process; }
0159   /// PDG particle number that identifies the type.
0160   constexpr Acts::PdgParticle pdg() const { return m_pdg; }
0161   /// Absolute PDG particle number that identifies the type.
0162   constexpr Acts::PdgParticle absolutePdg() const {
0163     return Acts::makeAbsolutePdgParticle(pdg());
0164   }
0165   /// Particle charge.
0166   constexpr Scalar charge() const { return m_charge; }
0167   /// Particle absolute charge.
0168   constexpr Scalar absoluteCharge() const { return std::abs(m_charge); }
0169   /// Particle mass.
0170   constexpr Scalar mass() const { return m_mass; }
0171
0172   /// Particle hypothesis.
0173   constexpr Acts::ParticleHypothesis hypothesis() const {
0174     return Acts::ParticleHypothesis(absolutePdg(), mass(), absoluteCharge());
0175   }
0176   /// Particl qOverP.
0177   constexpr Scalar qOverP() const {
0178     return hypothesis().qOverP(absoluteMomentum(), charge());
0179   }
0180
0181   /// Space-time position four-vector.
0182   constexpr const Vector4 &fourPosition() const { return m_position4; }
0183   /// Three-position, i.e. spatial coordinates without the time.
0184   auto position() const { return m_position4.segment<3>(Acts::ePos0); }
0185   /// Time coordinate.
0186   Scalar time() const { return m_position4[Acts::eTime]; }
0187   /// Energy-momentum four-vector.
0188   Vector4 fourMomentum() const {
0189     Vector4 mom4;
0190     // stored direction is always normalized
0191     mom4[Acts::eMom0] = m_absMomentum * m_direction[Acts::ePos0];
0192     mom4[Acts::eMom1] = m_absMomentum * m_direction[Acts::ePos1];
0193     mom4[Acts::eMom2] = m_absMomentum * m_direction[Acts::ePos2];
0194     mom4[Acts::eEnergy] = energy();
0195     return mom4;
0196   }
0197   /// Unit three-direction, i.e. the normalized momentum three-vector.
0198   const Vector3 &direction() const { return m_direction; }
0199   /// Polar angle.
0200   Scalar theta() const { return Acts::VectorHelpers::theta(direction()); }
0201   /// Azimuthal angle.
0202   Scalar phi() const { return Acts::VectorHelpers::phi(direction()); }
0203   /// Absolute momentum in the x-y plane.
0204   Scalar transverseMomentum() const {
0205     return m_absMomentum * m_direction.segment<2>(Acts::eMom0).norm();
0206   }
0207   /// Absolute momentum.
0208   constexpr Scalar absoluteMomentum() const { return m_absMomentum; }
0209   /// Absolute momentum.
0210   Vector3 momentum() const { return absoluteMomentum() * direction(); }
0211   /// Total energy, i.e. norm of the four-momentum.
0212   Scalar energy() const { return std::hypot(m_mass, m_absMomentum); }
0213
0214   /// Check if the particle is alive, i.e. is not at rest.
0215   constexpr bool isAlive() const { return Scalar{0} < m_absMomentum; }
0216
0217   constexpr bool isSecondary() const {
0218     return particleId().vertexSecondary() != 0 ||
0219            particleId().generation() != 0 || particleId().subParticle() != 0;
0220   }
0221
0222   // simulation specific properties
0223
0224   /// Set the proper time in the particle rest frame.
0225   ///
0226   /// @param properTime passed proper time in the rest frame
0227   constexpr Particle &setProperTime(Scalar properTime) {
0228     m_properTime = properTime;
0229     return *this;
0230   }
0231   /// Proper time in the particle rest frame.
0232   constexpr Scalar properTime() const { return m_properTime; }
0233
0234   /// Set the accumulated material measured in radiation/interaction lengths.
0235   ///
0236   /// @param pathInX0 accumulated material measured in radiation lengths
0237   /// @param pathInL0 accumulated material measured in interaction lengths
0238   constexpr Particle &setMaterialPassed(Scalar pathInX0, Scalar pathInL0) {
0239     m_pathInX0 = pathInX0;
0240     m_pathInL0 = pathInL0;
0241     return *this;
0242   }
0243   /// Accumulated path within material measured in radiation lengths.
0244   constexpr Scalar pathInX0() const { return m_pathInX0; }
0245   /// Accumulated path within material measured in interaction lengths.
0246   constexpr Scalar pathInL0() const { return m_pathInL0; }
0247
0248   /// Set the reference surface.
0249   ///
0250   /// @param surface reference surface
0251   Particle &setReferenceSurface(const Acts::Surface *surface) {
0252     m_referenceSurface = surface;
0253     return *this;
0254   }
0255
0256   /// Reference surface.
0257   const Acts::Surface *referenceSurface() const { return m_referenceSurface; }
0258
0259   /// Check if the particle has a reference surface.
0260   bool hasReferenceSurface() const { return m_referenceSurface != nullptr; }
0261
0262   /// Bound track parameters.
0263   Acts::Result<Acts::BoundTrackParameters> boundParameters(
0264       const Acts::GeometryContext &gctx) const {
0265     if (!hasReferenceSurface()) {
0266       return Acts::Result<Acts::BoundTrackParameters>::failure(
0267           std::error_code());
0268     }
0269     Acts::Result<Acts::Vector2> localResult =
0270         m_referenceSurface->globalToLocal(gctx, position(), direction());
0271     if (!localResult.ok()) {
0272       return localResult.error();
0273     }
0274     Acts::BoundVector params;
0275     params << localResult.value(), phi(), theta(), qOverP(), time();
0276     return Acts::BoundTrackParameters(referenceSurface()->getSharedPtr(),
0277                                       params, std::nullopt, hypothesis());
0278   }
0279
0280   Acts::CurvilinearTrackParameters curvilinearParameters() const {
0281     return Acts::CurvilinearTrackParameters(
0282         fourPosition(), direction(), qOverP(), std::nullopt, hypothesis());
0283   }
0284
0285   /// Set the number of hits.
0286   ///
0287   /// @param nHits number of hits
0288   constexpr Particle &setNumberOfHits(std::uint32_t nHits) {
0289     m_numberOfHits = nHits;
0290     return *this;
0291   }
0292
0293   /// Number of hits.
0294   constexpr std::uint32_t numberOfHits() const { return m_numberOfHits; }
0295
0296   /// Set the outcome of particle.
0297   ///
0298   /// @param outcome outcome code
0299   constexpr Particle &setOutcome(ParticleOutcome outcome) {
0300     m_outcome = outcome;
0301     return *this;
0302   }
0303
0304   /// Particle outcome.
0305   constexpr ParticleOutcome outcome() const { return m_outcome; }
0306
0307  private:
0308   // identity, i.e. things that do not change over the particle lifetime.
0309   /// Particle identifier within the event.
0310   Barcode m_particleId;
0311   /// Process type specifier.
0312   ProcessType m_process = ProcessType::eUndefined;
0313   /// PDG particle number.
0314   Acts::PdgParticle m_pdg = Acts::PdgParticle::eInvalid;
0315   // Particle charge and mass.
0316   Scalar m_charge = Scalar{0};
0317   Scalar m_mass = Scalar{0};
0318   // kinematics, i.e. things that change over the particle lifetime.
0319   Vector3 m_direction = Vector3::UnitZ();
0320   Scalar m_absMomentum = Scalar{0};
0321   Vector4 m_position4 = Vector4::Zero();
0322   /// proper time in the particle rest frame
0323   Scalar m_properTime = Scalar{0};
0324   // accumulated material
0325   Scalar m_pathInX0 = Scalar{0};
0326   Scalar m_pathInL0 = Scalar{0};
0327   /// number of hits
0328   std::uint32_t m_numberOfHits = 0;
0329   /// reference surface
0330   const Acts::Surface *m_referenceSurface{nullptr};
0331   /// outcome
0332   ParticleOutcome m_outcome = ParticleOutcome::Alive;
0333 };
0334
0335 std::ostream &operator<<(std::ostream &os, const Particle &particle);
0336
0337 }  // namespace ActsFatras