ActsFatras/EventData/Particle.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 "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/BoundTrackParameters.hpp"
0016 #include "Acts/EventData/ParticleHypothesis.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Utilities/MathHelpers.hpp"
0020 #include "Acts/Utilities/VectorHelpers.hpp"
0021 #include "ActsFatras/EventData/Barcode.hpp"
0022 #include "ActsFatras/EventData/GenerationProcess.hpp"
0023 #include "ActsFatras/EventData/SimulationOutcome.hpp"
0024
0025 #include <cmath>
0026 #include <iosfwd>
0027 #include <optional>
0028
0029 namespace ActsFatras {
0030
0031 /// Particle identity information and kinematic state.
0032 ///
0033 /// Also stores some simulation-specific properties.
0034 class Particle {
0035  public:
0036   /// Construct a default particle with invalid identity.
0037   Particle() = default;
0038   /// Construct a particle at rest with explicit mass and charge.
0039   ///
0040   /// @param particleId Particle identifier within an event
0041   /// @param pdg PDG id
0042   /// @param charge Particle charge in native units
0043   /// @param mass Particle mass in native units
0044   ///
0045   /// @warning It is the users responsibility that charge and mass match
0046   ///          the PDG particle number.
0047   Particle(Barcode particleId, Acts::PdgParticle pdg, double charge,
0048            double mass)
0049       : m_particleId(particleId), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
0050   /// Construct a particle at rest from a PDG particle number.
0051   ///
0052   /// @param particleId Particle identifier within an event
0053   /// @param pdg PDG particle number
0054   ///
0055   /// Charge and mass are retrieved from the particle data table.
0056   Particle(Barcode particleId, Acts::PdgParticle pdg);
0057   /// Copy constructor
0058   Particle(const Particle &) = default;
0059   /// Move constructor
0060   Particle(Particle &&) = default;
0061   /// Copy assignment operator
0062   /// @return Reference to this particle after copying
0063   Particle &operator=(const Particle &) = default;
0064   /// Move assignment operator
0065   /// @return Reference to this particle after moving
0066   Particle &operator=(Particle &&) = default;
0067
0068   /// Construct a new particle with a new identifier but same kinematics.
0069   ///
0070   /// @param particleId New particle identifier to assign
0071   /// @return New particle instance with updated identifier but same kinematics
0072   /// @note This is intentionally not a regular setter. The particle id
0073   ///       is used to identify the whole particle. Setting it on an existing
0074   ///       particle is usually a mistake.
0075   Particle withParticleId(Barcode particleId) const {
0076     Particle p = *this;
0077     p.m_particleId = particleId;
0078     return p;
0079   }
0080
0081   /// Kill the particle by setting the outcome to a non-alive status.
0082   /// @param outcome The outcome status to set for the killed particle (must not be Alive)
0083   /// @throws std::invalid_argument if the provided outcome is Alive
0084   void killParticle(SimulationOutcome outcome) {
0085     if (outcome == SimulationOutcome::Alive) {
0086       throw std::invalid_argument("Cannot kill particle with outcome 'Alive'.");
0087     }
0088     m_outcome = outcome;
0089   }
0090
0091   /// Set the process type that generated this particle.
0092   /// @param proc Process type that generated this particle
0093   /// @return Reference to this particle for method chaining
0094   Particle &setProcess(GenerationProcess proc) {
0095     m_process = proc;
0096     return *this;
0097   }
0098   /// Set the pdg.
0099   /// @param pdg PDG particle identifier
0100   /// @return Particle instance with updated PDG identifier
0101   Particle setPdg(Acts::PdgParticle pdg) {
0102     m_pdg = pdg;
0103     return *this;
0104   }
0105   /// Set the charge.
0106   /// @param charge Particle charge in native units
0107   /// @return Particle instance with updated charge
0108   Particle setCharge(double charge) {
0109     m_charge = charge;
0110     return *this;
0111   }
0112   /// Set the mass.
0113   /// @param mass Particle mass in native units
0114   /// @return Particle instance with updated mass
0115   Particle setMass(double mass) {
0116     m_mass = mass;
0117     return *this;
0118   }
0119   /// Set the particle ID.
0120   /// @param barcode New particle identifier barcode
0121   /// @return Reference to this particle for method chaining
0122   Particle &setParticleId(Barcode barcode) {
0123     m_particleId = barcode;
0124     return *this;
0125   }
0126   /// Set the space-time position four-vector.
0127   /// @param pos4 Four-vector containing spatial position and time
0128   /// @return Reference to this particle for method chaining
0129   Particle &setPosition4(const Acts::Vector4 &pos4) {
0130     m_position4 = pos4;
0131     return *this;
0132   }
0133   /// Set the space-time position four-vector from three-position and time.
0134   /// @param position Three-dimensional spatial position vector
0135   /// @param time Time coordinate
0136   /// @return Reference to this particle for method chaining
0137   Particle &setPosition4(const Acts::Vector3 &position, double time) {
0138     m_position4.segment<3>(Acts::ePos0) = position;
0139     m_position4[Acts::eTime] = time;
0140     return *this;
0141   }
0142   /// Set the space-time position four-vector from scalar components.
0143   /// @param x X coordinate
0144   /// @param y Y coordinate
0145   /// @param z Z coordinate
0146   /// @param time Time coordinate
0147   /// @return Reference to this particle for method chaining
0148   Particle &setPosition4(double x, double y, double z, double time) {
0149     m_position4[Acts::ePos0] = x;
0150     m_position4[Acts::ePos1] = y;
0151     m_position4[Acts::ePos2] = z;
0152     m_position4[Acts::eTime] = time;
0153     return *this;
0154   }
0155   /// Set the direction three-vector
0156   /// @param direction Three-dimensional direction vector (will be normalized)
0157   /// @return Reference to this particle for method chaining
0158   Particle &setDirection(const Acts::Vector3 &direction) {
0159     m_direction = direction;
0160     m_direction.normalize();
0161     return *this;
0162   }
0163   /// Set the direction three-vector from scalar components.
0164   /// @param dx X component of direction
0165   /// @param dy Y component of direction
0166   /// @param dz Z component of direction
0167   /// @return Reference to this particle for method chaining
0168   Particle &setDirection(double dx, double dy, double dz) {
0169     m_direction[Acts::ePos0] = dx;
0170     m_direction[Acts::ePos1] = dy;
0171     m_direction[Acts::ePos2] = dz;
0172     m_direction.normalize();
0173     return *this;
0174   }
0175   /// Set the absolute momentum.
0176   /// @param absMomentum Absolute momentum magnitude
0177   /// @return Reference to this particle for method chaining
0178   /// @throws std::invalid_argument if absMomentum is negative
0179   Particle &setAbsoluteMomentum(double absMomentum) {
0180     if (absMomentum < 0) {
0181       throw std::invalid_argument("Absolute momentum cannot be negative.");
0182     }
0183     m_absMomentum = absMomentum;
0184     return *this;
0185   }
0186
0187   /// Reduce the energy by the given amount. If the energy loss exceeds the
0188   /// current energy, the particle is killed with the specified outcome. If
0189   /// stopping the particle was not expected (stoppedOutcome is Alive), an
0190   /// exception is thrown.
0191   /// @param delta The energy loss amount to subtract from the current energy
0192   /// @param stoppedOutcome The outcome to set if the energy loss exceeds the current energy
0193   /// @return Reference to this particle for method chaining
0194   /// @throws std::invalid_argument if the energy loss exceeds the current energy and stoppedOutcome is Alive
0195   Particle &loseEnergy(double delta, SimulationOutcome stoppedOutcome =
0196                                          SimulationOutcome::Alive) {
0197     const double newEnergy = energy() - delta;
0198     if (newEnergy <= m_mass) {
0199       if (stoppedOutcome == SimulationOutcome::Alive) {
0200         throw std::invalid_argument(
0201             "Energy loss cannot exceed the current energy of the particle if "
0202             "the particle is to remain alive.");
0203       }
0204       killParticle(stoppedOutcome);
0205     } else {
0206       m_absMomentum = Acts::fastCathetus(newEnergy, m_mass);
0207     }
0208     return *this;
0209   }
0210
0211   /// Particle identifier within an event.
0212   /// @return The unique particle identifier barcode
0213   Barcode particleId() const { return m_particleId; }
0214   /// Which type of process generated this particle.
0215   /// @return The process type that generated this particle
0216   GenerationProcess process() const { return m_process; }
0217   /// PDG particle number that identifies the type.
0218   /// @return The PDG particle identifier
0219   Acts::PdgParticle pdg() const { return m_pdg; }
0220   /// Absolute PDG particle number that identifies the type.
0221   /// @return The absolute PDG particle identifier (positive value)
0222   Acts::PdgParticle absolutePdg() const {
0223     return Acts::makeAbsolutePdgParticle(pdg());
0224   }
0225   /// Particle charge.
0226   /// @return The particle charge in native units
0227   double charge() const { return m_charge; }
0228   /// Particle absolute charge.
0229   /// @return The absolute particle charge (positive value)
0230   double absoluteCharge() const { return std::abs(m_charge); }
0231   /// Particle mass.
0232   /// @return The particle mass in native units
0233   double mass() const { return m_mass; }
0234
0235   /// Particle hypothesis.
0236   /// @return Particle hypothesis containing PDG, mass, and charge information
0237   Acts::ParticleHypothesis hypothesis() const {
0238     return Acts::ParticleHypothesis(
0239         absolutePdg(), static_cast<float>(mass()),
0240         Acts::ChargeHypothesis{static_cast<float>(absoluteCharge())},
0241         charge() != 0 ? std::nullopt
0242                       : std::optional<double>(absoluteMomentum()));
0243   }
0244   /// Particl qOverP.
0245   /// @return The charge over momentum ratio
0246   double qOverP() const {
0247     return hypothesis().qOverP(absoluteMomentum(), charge());
0248   }
0249
0250   /// Space-time position four-vector.
0251   /// @return Reference to the four-dimensional position vector (x, y, z, t)
0252   const Acts::Vector4 &fourPosition() const { return m_position4; }
0253   /// Three-position, i.e. spatial coordinates without the time.
0254   /// @return Three-dimensional position vector (x, y, z)
0255   auto position() const { return m_position4.segment<3>(Acts::ePos0); }
0256   /// Time coordinate.
0257   /// @return The time coordinate value
0258   double time() const { return m_position4[Acts::eTime]; }
0259   /// Energy-momentum four-vector.
0260   /// @return Four-dimensional momentum vector (px, py, pz, E)
0261   Acts::Vector4 fourMomentum() const {
0262     Acts::Vector4 mom4;
0263     // stored direction is always normalized
0264     mom4[Acts::eMom0] = m_absMomentum * m_direction[Acts::ePos0];
0265     mom4[Acts::eMom1] = m_absMomentum * m_direction[Acts::ePos1];
0266     mom4[Acts::eMom2] = m_absMomentum * m_direction[Acts::ePos2];
0267     mom4[Acts::eEnergy] = energy();
0268     return mom4;
0269   }
0270   /// Unit three-direction, i.e. the normalized momentum three-vector.
0271   /// @return Reference to the normalized direction vector
0272   const Acts::Vector3 &direction() const { return m_direction; }
0273   /// Polar angle.
0274   /// @return The polar angle (theta) in radians
0275   double theta() const { return Acts::VectorHelpers::theta(direction()); }
0276   /// Azimuthal angle.
0277   /// @return The azimuthal angle (phi) in radians
0278   double phi() const { return Acts::VectorHelpers::phi(direction()); }
0279   /// Absolute momentum in the x-y plane.
0280   /// @return The transverse momentum magnitude
0281   double transverseMomentum() const {
0282     return m_absMomentum * m_direction.segment<2>(Acts::eMom0).norm();
0283   }
0284   /// Absolute momentum.
0285   /// @return The absolute momentum magnitude
0286   double absoluteMomentum() const { return m_absMomentum; }
0287   /// Absolute momentum.
0288   /// @return Three-dimensional momentum vector
0289   Acts::Vector3 momentum() const { return absoluteMomentum() * direction(); }
0290   /// Total energy, i.e. norm of the four-momentum.
0291   /// @return The total energy calculated from mass and momentum
0292   double energy() const { return std::hypot(m_mass, m_absMomentum); }
0293
0294   /// Check if the particle is alive, i.e. is not at rest.
0295   /// @return True if particle has non-zero momentum, false otherwise
0296   bool isAlive() const { return m_outcome == SimulationOutcome::Alive; }
0297
0298   /// Check if this is a secondary particle.
0299   /// @return True if particle is a secondary (has non-zero vertex secondary, generation, or sub-particle), false otherwise
0300   bool isSecondary() const {
0301     return particleId().vertexSecondary() != 0 ||
0302            particleId().generation() != 0 || particleId().subParticle() != 0;
0303   }
0304
0305   // simulation specific properties
0306
0307   /// Set the proper time in the particle rest frame.
0308   ///
0309   /// @param properTime passed proper time in the rest frame
0310   /// @return Reference to this particle for method chaining
0311   Particle &setProperTime(double properTime) {
0312     m_properTime = properTime;
0313     return *this;
0314   }
0315   /// Proper time in the particle rest frame.
0316   /// @return The proper time in the rest frame
0317   double properTime() const { return m_properTime; }
0318
0319   /// Set the accumulated material measured in radiation/interaction lengths.
0320   ///
0321   /// @param pathInX0 accumulated material measured in radiation lengths
0322   /// @param pathInL0 accumulated material measured in interaction lengths
0323   /// @return Reference to this particle for method chaining
0324   Particle &setMaterialPassed(double pathInX0, double pathInL0) {
0325     m_pathInX0 = pathInX0;
0326     m_pathInL0 = pathInL0;
0327     return *this;
0328   }
0329   /// Accumulated path within material measured in radiation lengths.
0330   /// @return The accumulated path in radiation lengths
0331   double pathInX0() const { return m_pathInX0; }
0332   /// Accumulated path within material measured in interaction lengths.
0333   /// @return The accumulated path in interaction lengths
0334   double pathInL0() const { return m_pathInL0; }
0335
0336   /// Set the reference surface.
0337   ///
0338   /// @param surface Reference surface for bound track parameters
0339   /// @return Reference to this particle for method chaining
0340   Particle &setReferenceSurface(const Acts::Surface *surface) {
0341     m_referenceSurface = surface;
0342     return *this;
0343   }
0344
0345   /// Reference surface.
0346   /// @return Pointer to the reference surface, or nullptr if not set
0347   const Acts::Surface *referenceSurface() const { return m_referenceSurface; }
0348
0349   /// Check if the particle has a reference surface.
0350   /// @return True if reference surface is set, false otherwise
0351   bool hasReferenceSurface() const { return m_referenceSurface != nullptr; }
0352
0353   /// Bound track parameters.
0354   /// @param gctx Geometry context for coordinate transformations
0355   /// @return Result containing bound track parameters or error if no reference surface
0356   Acts::Result<Acts::BoundTrackParameters> boundParameters(
0357       const Acts::GeometryContext &gctx) const {
0358     if (!hasReferenceSurface()) {
0359       return Acts::Result<Acts::BoundTrackParameters>::failure(
0360           std::error_code());
0361     }
0362     Acts::Result<Acts::Vector2> localResult =
0363         m_referenceSurface->globalToLocal(gctx, position(), direction());
0364     if (!localResult.ok()) {
0365       return localResult.error();
0366     }
0367     Acts::BoundVector params;
0368     params << localResult.value(), phi(), theta(), qOverP(), time();
0369     return Acts::BoundTrackParameters(referenceSurface()->getSharedPtr(),
0370                                       params, std::nullopt, hypothesis());
0371   }
0372
0373   /// @return Curvilinear track parameters representation
0374   Acts::BoundTrackParameters curvilinearParameters() const {
0375     return Acts::BoundTrackParameters::createCurvilinear(
0376         fourPosition(), direction(), qOverP(), std::nullopt, hypothesis());
0377   }
0378
0379   /// Set the number of hits.
0380   ///
0381   /// @param nHits number of hits
0382   /// @return Reference to this particle for method chaining
0383   Particle &setNumberOfHits(std::uint32_t nHits) {
0384     m_numberOfHits = nHits;
0385     return *this;
0386   }
0387
0388   /// Number of hits.
0389   /// @return The number of hits associated with this particle
0390   std::uint32_t numberOfHits() const { return m_numberOfHits; }
0391
0392   /// Set the parent particle id.
0393   ///
0394   /// A default-constructed @c Barcode means "unknown parent". Populated by
0395   /// input converters that have access to parent relationships (e.g. EDM4hep);
0396   /// left at the default by simulation engines that don't track it.
0397   ///
0398   /// @param parentId Barcode identifying the parent particle
0399   /// @return Reference to this particle for method chaining
0400   Particle &setParentParticleId(Barcode parentId) {
0401     m_parentParticleId = parentId;
0402     return *this;
0403   }
0404
0405   /// Parent particle id, or default-constructed @c Barcode if unknown.
0406   /// @return The parent particle barcode
0407   Barcode parentParticleId() const { return m_parentParticleId; }
0408
0409   /// Set the outcome of particle.
0410   ///
0411   /// @param outcome outcome code
0412   /// @return Reference to this particle for method chaining
0413   Particle &setOutcome(SimulationOutcome outcome) {
0414     m_outcome = outcome;
0415     return *this;
0416   }
0417
0418   /// Particle outcome.
0419   /// @return The outcome status of this particle
0420   SimulationOutcome outcome() const { return m_outcome; }
0421
0422  private:
0423   // identity, i.e. things that do not change over the particle lifetime.
0424   /// Particle identifier within the event.
0425   Barcode m_particleId;
0426   /// Process type specifier.
0427   GenerationProcess m_process = GenerationProcess::eUndefined;
0428   /// PDG particle number.
0429   Acts::PdgParticle m_pdg = Acts::PdgParticle::eInvalid;
0430   // Particle charge and mass.
0431   double m_charge = 0;
0432   double m_mass = 0;
0433   // kinematics, i.e. things that change over the particle lifetime.
0434   Acts::Vector3 m_direction = Acts::Vector3::UnitZ();
0435   double m_absMomentum = 0;
0436   Acts::Vector4 m_position4 = Acts::Vector4::Zero();
0437   /// proper time in the particle rest frame
0438   double m_properTime = 0;
0439   // accumulated material
0440   double m_pathInX0 = 0;
0441   double m_pathInL0 = 0;
0442   /// number of hits
0443   std::uint32_t m_numberOfHits = 0;
0444   /// parent particle id (default-constructed = unknown)
0445   Barcode m_parentParticleId{};
0446   /// reference surface
0447   const Acts::Surface *m_referenceSurface{nullptr};
0448   /// outcome
0449   SimulationOutcome m_outcome = SimulationOutcome::Alive;
0450 };
0451
0452 /// Print particle to output stream
0453 /// @param os Output stream
0454 /// @param particle Particle to print
0455 /// @return Output stream
0456 std::ostream &operator<<(std::ostream &os, const Particle &particle);
0457
0458 }  // namespace ActsFatras