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-2020 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/Utilities/Definitions.hpp"
#include "Acts/Utilities/PdgParticle.hpp"
#include "ActsFatras/EventData/Barcode.hpp"
#include "ActsFatras/EventData/ProcessType.hpp"

#include <cmath>
#include <iosfwd>
#include <limits>

namespace ActsFatras {

class Particle {
 public:
  using Scalar = double;
  using Vector3 = Acts::ActsVector<Scalar, 3>;
  using Vector4 = Acts::ActsVector<Scalar, 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) { return m_process = proc, *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_unitDirection = direction;
    m_unitDirection.normalize();
    return *this;
  }
  Particle &setDirection(Scalar dx, Scalar dy, Scalar dz) {
    m_unitDirection[Acts::ePos0] = dx;
    m_unitDirection[Acts::ePos1] = dy;
    m_unitDirection[Acts::ePos2] = dz;
    m_unitDirection.normalize();
    return *this;
  }
  Particle &setAbsMomentum(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 Scalar charge() const { return m_charge; }
  constexpr Scalar mass() const { return m_mass; }

  constexpr const Vector4 &position4() const { return m_position4; }
  auto position() const { return m_position4.segment<3>(Acts::ePos0); }
  Scalar time() const { return m_position4[Acts::eTime]; }
  Vector4 momentum4() const {
    Vector4 mom4;
    // stored direction is always normalized
    mom4[Acts::eMom0] = m_absMomentum * m_unitDirection[Acts::ePos0];
    mom4[Acts::eMom1] = m_absMomentum * m_unitDirection[Acts::ePos1];
    mom4[Acts::eMom2] = m_absMomentum * m_unitDirection[Acts::ePos2];
    mom4[Acts::eEnergy] = energy();
    return mom4;
  }
  const Vector3 &unitDirection() const { return m_unitDirection; }
  Scalar transverseMomentum() const {
    return m_absMomentum * m_unitDirection.segment<2>(Acts::eMom0).norm();
  }
  constexpr Scalar absMomentum() const { return m_absMomentum; }
  Scalar energy() const { return std::hypot(m_mass, m_absMomentum); }

  constexpr operator bool() const { return Scalar(0) < m_absMomentum; }
  constexpr bool operator!() const { return m_absMomentum <= Scalar(0); }

  constexpr Particle &setMaterialPassed(Scalar pathX0, Scalar pathL0) {
    m_pathX0 = pathX0;
    m_pathL0 = pathL0;
    return *this;
  }
  constexpr Particle &setMaterialLimits(Scalar limitX0, Scalar limitL0) {
    m_limitX0 = limitX0;
    m_limitL0 = limitL0;
    return *this;
  }
  constexpr Scalar pathInX0() const { return m_pathX0; }
  constexpr Scalar pathInL0() const { return m_pathL0; }
  constexpr Scalar pathLimitX0() const { return m_limitX0; }
  constexpr Scalar pathLimitL0() const { return m_limitL0; }

 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_unitDirection = Vector3::UnitZ();
  Scalar m_absMomentum = Scalar(0);
  Vector4 m_position4 = Vector4::Zero();
  // simulation-specific X0/L0 information and limits
  // these values are here to simplify the simulation of (nuclear) interactions.
  // instead of checking at every surface whether an interaction should occur we
  // can draw an overall limit once. the relevant interaction only needs to
  // be executed once the limit is reached.
  // this information is not really particle-specific and should probably be
  // handled separately. for now, storing it directly here is the simplest
  // solution.
  Scalar m_pathX0 = Scalar(0);
  Scalar m_pathL0 = Scalar(0);
  Scalar m_limitX0 = std::numeric_limits<Scalar>::max();
  Scalar m_limitL0 = std::numeric_limits<Scalar>::max();
};

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

}  // namespace ActsFatras