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