Particle.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Particle.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2018-2021 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Definitions/PdgParticle.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/EventData/ParticleHypothesis.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/VectorHelpers.hpp"
#include "ActsFatras/EventData/Barcode.hpp"
#include "ActsFatras/EventData/ProcessType.hpp"

#include <algorithm>
#include <cmath>
#include <iosfwd>
#include <limits>
#include <memory>
#include <optional>

namespace ActsFatras {

class Particle {
 public:
  using Scalar = Acts::ActsScalar;
  using Vector3 = Acts::ActsVector<3>;
  using Vector4 = Acts::ActsVector<4>;

  Particle() = default;
  Particle(Barcode particleId, Acts::PdgParticle pdg, Scalar charge,
           Scalar mass)
      : m_particleId(particleId), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
  Particle(Barcode particleId, Acts::PdgParticle pdg);
  Particle(const Particle &) = default;
  Particle(Particle &&) = default;
  Particle &operator=(const Particle &) = default;
  Particle &operator=(Particle &&) = default;

  Particle withParticleId(Barcode particleId) const {
    Particle p = *this;
    p.m_particleId = particleId;
    return p;
  }

  Particle &setProcess(ProcessType proc) {
    m_process = proc;
    return *this;
  }
  Particle setPdg(Acts::PdgParticle pdg) {
    m_pdg = pdg;
    return *this;
  }
  Particle setCharge(Scalar charge) {
    m_charge = charge;
    return *this;
  }
  Particle setMass(Scalar mass) {
    m_mass = mass;
    return *this;
  }
  Particle &setParticleId(Barcode barcode) {
    m_particleId = barcode;
    return *this;
  }
  Particle &setPosition4(const Vector4 &pos4) {
    m_position4 = pos4;
    return *this;
  }
  Particle &setPosition4(const Vector3 &position, Scalar time) {
    m_position4.segment<3>(Acts::ePos0) = position;
    m_position4[Acts::eTime] = time;
    return *this;
  }
  Particle &setPosition4(Scalar x, Scalar y, Scalar z, Scalar time) {
    m_position4[Acts::ePos0] = x;
    m_position4[Acts::ePos1] = y;
    m_position4[Acts::ePos2] = z;
    m_position4[Acts::eTime] = time;
    return *this;
  }
  Particle &setDirection(const Vector3 &direction) {
    m_direction = direction;
    m_direction.normalize();
    return *this;
  }
  Particle &setDirection(Scalar dx, Scalar dy, Scalar dz) {
    m_direction[Acts::ePos0] = dx;
    m_direction[Acts::ePos1] = dy;
    m_direction[Acts::ePos2] = dz;
    m_direction.normalize();
    return *this;
  }
  Particle &setAbsoluteMomentum(Scalar absMomentum) {
    m_absMomentum = absMomentum;
    return *this;
  }

  Particle &correctEnergy(Scalar delta) {
    const auto newEnergy = std::hypot(m_mass, m_absMomentum) + delta;
    if (newEnergy <= m_mass) {
      m_absMomentum = Scalar(0);
    } else {
      m_absMomentum = std::sqrt(newEnergy * newEnergy - m_mass * m_mass);
    }
    return *this;
  }

  constexpr Barcode particleId() const { return m_particleId; }
  constexpr ProcessType process() const { return m_process; }
  constexpr Acts::PdgParticle pdg() const { return m_pdg; }
  constexpr Acts::PdgParticle absolutePdg() const {
    return Acts::makeAbsolutePdgParticle(pdg());
  }
  constexpr Scalar charge() const { return m_charge; }
  constexpr Scalar absoluteCharge() const { return std::abs(m_charge); }
  constexpr Scalar mass() const { return m_mass; }

  constexpr Acts::ParticleHypothesis hypothesis() const {
    return Acts::ParticleHypothesis(absolutePdg(), mass(), absoluteCharge());
  }
  constexpr Scalar qOverP() const {
    return hypothesis().qOverP(absoluteMomentum(), charge());
  }

  constexpr const Vector4 &fourPosition() const { return m_position4; }
  auto position() const { return m_position4.segment<3>(Acts::ePos0); }
  Scalar time() const { return m_position4[Acts::eTime]; }
  Vector4 fourMomentum() const {
    Vector4 mom4;
    // stored direction is always normalized
    mom4[Acts::eMom0] = m_absMomentum * m_direction[Acts::ePos0];
    mom4[Acts::eMom1] = m_absMomentum * m_direction[Acts::ePos1];
    mom4[Acts::eMom2] = m_absMomentum * m_direction[Acts::ePos2];
    mom4[Acts::eEnergy] = energy();
    return mom4;
  }
  const Vector3 &direction() const { return m_direction; }
  Scalar theta() const { return Acts::VectorHelpers::theta(direction()); }
  Scalar phi() const { return Acts::VectorHelpers::phi(direction()); }
  Scalar transverseMomentum() const {
    return m_absMomentum * m_direction.segment<2>(Acts::eMom0).norm();
  }
  constexpr Scalar absoluteMomentum() const { return m_absMomentum; }
  Vector3 momentum() const { return absoluteMomentum() * direction(); }
  Scalar energy() const { return std::hypot(m_mass, m_absMomentum); }

  constexpr bool isAlive() const { return Scalar(0) < m_absMomentum; }

  constexpr bool isSecondary() const {
    return particleId().vertexSecondary() != 0 ||
           particleId().generation() != 0 || particleId().subParticle() != 0;
  }

  // simulation specific properties

  constexpr Particle &setProperTime(Scalar properTime) {
    m_properTime = properTime;
    return *this;
  }
  constexpr Scalar properTime() const { return m_properTime; }

  constexpr Particle &setMaterialPassed(Scalar pathInX0, Scalar pathInL0) {
    m_pathInX0 = pathInX0;
    m_pathInL0 = pathInL0;
    return *this;
  }
  constexpr Scalar pathInX0() const { return m_pathInX0; }
  constexpr Scalar pathInL0() const { return m_pathInL0; }

  Particle &setReferenceSurface(const Acts::Surface *surface) {
    m_referenceSurface = surface;
    return *this;
  }

  const Acts::Surface *referenceSurface() const { return m_referenceSurface; }

  bool hasReferenceSurface() const { return m_referenceSurface != nullptr; }

  Acts::Result<Acts::BoundTrackParameters> boundParameters(
      const Acts::GeometryContext &gctx) const {
    if (!hasReferenceSurface()) {
      return Acts::Result<Acts::BoundTrackParameters>::failure(
          std::error_code());
    }
    Acts::Result<Acts::Vector2> localResult =
        m_referenceSurface->globalToLocal(gctx, position(), direction());
    if (!localResult.ok()) {
      return localResult.error();
    }
    Acts::BoundVector params;
    params << localResult.value(), phi(), theta(), qOverP(), time();
    return Acts::BoundTrackParameters(referenceSurface()->getSharedPtr(),
                                      params, std::nullopt, hypothesis());
  }

  Acts::CurvilinearTrackParameters curvilinearParameters() const {
    return Acts::CurvilinearTrackParameters(
        fourPosition(), direction(), qOverP(), std::nullopt, hypothesis());
  }

 private:
  // identity, i.e. things that do not change over the particle lifetime.
  Barcode m_particleId;
  ProcessType m_process = ProcessType::eUndefined;
  Acts::PdgParticle m_pdg = Acts::PdgParticle::eInvalid;
  // Particle charge and mass.
  Scalar m_charge = Scalar(0);
  Scalar m_mass = Scalar(0);
  // kinematics, i.e. things that change over the particle lifetime.
  Vector3 m_direction = Vector3::UnitZ();
  Scalar m_absMomentum = Scalar(0);
  Vector4 m_position4 = Vector4::Zero();
  // proper time in the particle rest frame
  Scalar m_properTime = Scalar(0);
  // accumulated material
  Scalar m_pathInX0 = Scalar(0);
  Scalar m_pathInL0 = Scalar(0);
  // reference surface
  const Acts::Surface *m_referenceSurface{nullptr};
};

std::ostream &operator<<(std::ostream &os, const Particle &particle);

}  // namespace ActsFatras