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