ParticleContainer.hpp

Go to the documentation of this file.

 
#ifndef OPAL_PARTICLE_CONTAINER_H
#define OPAL_PARTICLE_CONTAINER_H
 
// #include <functional>
#include <cmath>
#include <cstddef>
#include <functional>
#include <memory>
#include <string>
#include <vector>
 
#include "Manager/BaseManager.h"
 
#include "PartBunch/FieldContainer.hpp"
 
#include "Algorithms/CoordinateSystemTrafo.h"
#include "Algorithms/DistributionMoments.h"
#include "Algorithms/PartData.h"
#include "Algorithms/Quaternion.hpp"
#include "PartBunch/BunchStateHandler.h"
 
#include "Utilities/OpalException.h"
#include "Utilities/Options.h"
 
#include "Physics/Physics.h"
 
// #include <Kokkos_Core.hpp>
 
template <typename T>
using ParticleAttrib = ippl::ParticleAttrib<T>;
 
using size_type = ippl::detail::size_type;
 
#include "Processes/GlobalProcesses/GlobalProcess.h"
 
template <typename T, unsigned Dim = 3>
class ParticleContainer
    : public ippl::ParticleBase<
              ippl::ParticleSpatialLayout<T, Dim>, Kokkos::DefaultExecutionSpace::memory_space> {
    using Base = ippl::ParticleBase<
            ippl::ParticleSpatialLayout<T, Dim>, Kokkos::DefaultExecutionSpace::memory_space>;
 
private:
    using Base::alloc;
    using Base::create;
    using Base::destroy;
 
public:
    enum class QMStorageMode { SingleValue, Attributes };
 
    using bin_index_type = short int;  // Needed in AdaptBins class
 
    using qm_view_type = typename ippl::ParticleAttrib<double>::view_type;
 
    using spin_vector_type = ippl::Vector<float, 3>;
 
public:
    qm_view_type getQView() const {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            return QAttr.getView();
        }
        return QView_m;
    }
 
    qm_view_type getMView() const {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            return MAttr.getView();
        }
        return MView_m;
    }
 
    ippl::ParticleAttrib<double> dt;
 
    ippl::ParticleAttrib<double> Phi;
 
    ippl::ParticleAttrib<bin_index_type> Bin;
 
    short Sp = 0;
 
    typename Base::particle_position_type P;
 
    typename Base::particle_position_type E;
 
    // typename Base::particle_position_type Etmp; // TODO: might not need this...
 
    typename Base::particle_position_type B;
 
    ippl::ParticleAttrib<bool> InvalidMask;
 
    ippl::ParticleAttrib<spin_vector_type> Pol;
 
    bool hasSpin() const { return spinEnabled_m; }
 
    ParticleContainer(Mesh_t<Dim>& mesh, FieldLayout_t<Dim>& FL, bool spinEnabled = false)
        : pl_m(FL, mesh),
          qmStorageMode_m(
                  Options::useQMAttributes ? QMStorageMode::Attributes
                                           : QMStorageMode::SingleValue),
          distMoments_m(),
          QView_m("ParticleContainer::QView_m", 1),
          MView_m("ParticleContainer::MView_m", 1),
          spinEnabled_m(spinEnabled) {
        this->initialize(pl_m);
        registerAttributes();
        setupBCs();
        Kokkos::deep_copy(QView_m, 0.0);
        Kokkos::deep_copy(MView_m, 0.0);
    }
 
    ~ParticleContainer() {}
 
    void registerAttributes() {
        // register the particle attributes
        this->addAttribute(dt);
        this->addAttribute(Phi);
        this->addAttribute(Bin);
        this->addAttribute(P);
        this->addAttribute(E);
        this->addAttribute(B);
        this->addAttribute(InvalidMask);
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            this->addAttribute(QAttr);
            this->addAttribute(MAttr);
        }
        if (spinEnabled_m) {
            this->addAttribute(Pol);
        }
    }
 
    void setupBCs() { setBCAllPeriodic(); }
 
    void transformBunch(const CoordinateSystemTrafo& trafo) {
        const size_t nLoc = this->getLocalNum();
        trafo.transformBunchTo(this->R.getView(), nLoc);
        trafo.rotateBunchTo(this->P.getView(), nLoc);
        trafo.rotateBunchTo(this->E.getView(), nLoc);
        trafo.rotateBunchTo(this->B.getView(), nLoc);
    }
    PLayout_t<T, Dim>& getPL() { return pl_m; }
 
    void setBunchStateHandler(std::shared_ptr<BunchStateHandler> handler) {
        // We only keep the slot: per-container flags own their own sync, so
        // no back-reference to the handler is needed. The handler itself
        // manages bunch-wide state, which ParticleContainer never touches.
        containerState_m = handler->registerContainer();
        distMoments_m.setContainerState(containerState_m);
    }
 
    // -- per-container state pass-throughs --------------------------------
    // Thin wrappers around the slot's own API so callers don't need to reach
    // into `containerState_m` directly.
 
    bool isUnitlessPositions() const { return containerState_m->unitlessPositions; }
 
    bool isMomentsDirty() const { return containerState_m->momentsDirty; }
    void markMomentsDirty() { containerState_m->markMomentsDirty(); }
    void markMomentsClean() { containerState_m->markMomentsClean(); }
 
    void updateMoments() {
        /*
        Quick check that this container was registered with a BunchStateHandler. If not, throw an
        error. This seems to be an easy to make error when interacting (e.g. through unit tests)
        with the ParticleContainer.
        */
        if (!containerState_m) {
            throw OpalException(
                    "ParticleContainer::updateMoments",
                    "BunchStateHandler not set in ParticleContainer (containerState is null).");
        }
 
        size_t Np = this->getTotalNum();
        Np = (Np == 0) ? 1 : Np;  // only used for normalization in the moments class --> avoid
                                  // division by zero
 
        size_t Nlocal = this->getLocalNum();
 
        distMoments_m.computeMoments(this->R.getView(), this->P.getView(), getMView(), Np, Nlocal);
    }
 
    void setEnergyReferenceMass(double referenceMassGeV, bool rescaleToReference = true) {
        distMoments_m.setEnergyReferenceMass(referenceMassGeV, rescaleToReference);
    }
 
    Vector_t<double, 3> getMeanP() const { return distMoments_m.getMeanMomentum(); }
 
    Vector_t<double, 3> getRmsP() const { return distMoments_m.getStandardDeviationMomentum(); }
 
    Vector_t<double, 3> getMeanR() const { return distMoments_m.getMeanPosition(); }
 
    Vector_t<double, 3> getRmsR() const { return distMoments_m.getStandardDeviationPosition(); }
 
    Vector_t<double, 3> getRmsRP() const { return distMoments_m.getStandardDeviationRP(); }
 
    void computeMinMaxR() {
        size_t Nlocal = this->getLocalNum();
        distMoments_m.computeMinMaxPosition(this->R.getView(), Nlocal);
    }
 
    Vector_t<double, 3> getMinR() const { return distMoments_m.getMinPosition(); }
 
    Vector_t<double, 3> getMaxR() const { return distMoments_m.getMaxPosition(); }
 
    matrix6x6_t getCovMatrix() const { return distMoments_m.getMoments6x6(); }
 
    double getMeanKineticEnergy() const { return distMoments_m.getMeanKineticEnergy(); }
 
    double getStdKineticEnergy() const { return distMoments_m.getStdKineticEnergy(); }
 
    Vector_t<double, 6> getMeans() const { return distMoments_m.getMeans(); }
 
    Vector_t<double, 6> getCentroid() const { return distMoments_m.getCentroid(); }
 
    Vector_t<double, 3> getNormEmit() const { return distMoments_m.getNormalizedEmittance(); }
 
    Vector_t<double, 3> getGeometricEmit() const { return distMoments_m.getGeometricEmittance(); }
 
    double getDx() const { return distMoments_m.getDx(); }
 
    double getDDx() const { return distMoments_m.getDDx(); }
 
    double getDy() const { return distMoments_m.getDy(); }
 
    double getDDy() const { return distMoments_m.getDDy(); }
 
    double getDebyeLength() const { return distMoments_m.getDebyeLength(); }
 
    double getMeanGammaZ() const { return distMoments_m.getMeanGammaZ(); }
 
    double getTemperature() const { return distMoments_m.getTemperature(); }
 
    double getPlasmaParameter() const { return distMoments_m.getPlasmaParameter(); }
 
    double computeDebyeLength(double density) {
        size_t Np = this->getTotalNum();
        Np = (Np == 0) ? 1 : Np;  // only used for normalization in the moments class --> avoid
                                  // division by zero
 
        size_t Nlocal = this->getLocalNum();
        distMoments_m.computeDebyeLength(this->P.getView(), Np, Nlocal, density);
        return distMoments_m.getDebyeLength();
    }
 
    void setQ(double q) {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto view         = QAttr.getView();
            const size_type n = this->getLocalNum();
            if (n == 0) {
                return;
            }
            Kokkos::parallel_for(
                    "ParticleContainer::setQ", n,
                    KOKKOS_LAMBDA(const size_type i) { view(i) = q; });
            Kokkos::fence();
        } else {
            Kokkos::deep_copy(QView_m, q);
        }
    }
 
    double getChargePerParticle() const {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto view = QAttr.getView();
            if (view.extent(0) == 0) {
                return 0.0;
            }
            double q = 0.0;
            Kokkos::deep_copy(q, Kokkos::subview(view, 0));
            return q;
        }
        double q = 0.0;
        Kokkos::deep_copy(q, Kokkos::subview(QView_m, 0));
        return q;
    }
 
    double getTotalCharge() const { return getChargePerParticle() * this->getTotalNum(); }
 
    void setM(double m) {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto view         = MAttr.getView();
            const size_type n = view.extent(0);
 
            Kokkos::parallel_for(
                    "ParticleContainer::setM", n,
                    KOKKOS_LAMBDA(const size_type i) { view(i) = m; });
            Kokkos::fence();
        } else {
            Kokkos::deep_copy(MView_m, m);
        }
    }
 
    double getMassPerParticle() const {
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto view = MAttr.getView();
            if (view.extent(0) == 0) {
                return 0.0;
            }
            double m = 0.0;
            Kokkos::deep_copy(m, Kokkos::subview(view, 0));
            return m;
        }
        double m = 0.0;
        Kokkos::deep_copy(m, Kokkos::subview(MView_m, 0));
        return m;
    }
 
    double getTotalMass() const { return getMassPerParticle() * this->getTotalNum(); }
 
    const Vector_t<double, Dim>& getRefPartR() const { return refPartR_m; }
 
    Vector_t<double, Dim>& getRefPartR() { return refPartR_m; }
 
    void setRefPartR(const Vector_t<double, Dim>& refPartR) { refPartR_m = refPartR; }
 
    const Vector_t<double, Dim>& getRefPartP() const { return refPartP_m; }
 
    Vector_t<double, Dim>& getRefPartP() { return refPartP_m; }
 
    void setRefPartP(const Vector_t<double, Dim>& refPartP) { refPartP_m = refPartP; }
 
    void setReference(const PartData* ref) {
        reference_m = ref;
        if (reference_m) {
            // PartData mass is stored in eV; DistributionMoments expects GeV.
            setEnergyReferenceMass(reference_m->getM() * Units::eV2GeV, true);
        }
    }
 
    const PartData* getReference() const { return reference_m; }
 
    void set_sPos(double sPos) { sPos_m = sPos; }
 
    double get_sPos() const { return sPos_m; }
 
    void setGlobalToLocalQuaternion(const Quaternion_t& globalToLocalQuaternion) {
        globalToLocalQuaternion_m = globalToLocalQuaternion;
    }
 
    Quaternion_t getGlobalToLocalQuaternion() const { return globalToLocalQuaternion_m; }
 
    const CoordinateSystemTrafo& getToLabTrafo() const { return toLabTrafo_m; }
 
    CoordinateSystemTrafo& getToLabTrafo() { return toLabTrafo_m; }
 
    void setToLabTrafo(const CoordinateSystemTrafo& toLabTrafo) { toLabTrafo_m = toLabTrafo; }
 
    void updateRefToLabCSTrafo(double bunchDT) {
        Vector_t<double, 3> R = toLabTrafo_m.transformFrom(refPartR_m);
        Vector_t<double, 3> P = toLabTrafo_m.rotateFrom(refPartP_m);
 
        const double ds =
                std::copysign(1.0, bunchDT) * std::sqrt(R[0] * R[0] + R[1] * R[1] + R[2] * R[2]);
        sPos_m += ds;
 
        CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t<double, 3>(0, 0, 1)));
        transformBunch(update);
        toLabTrafo_m = toLabTrafo_m * update.inverted();
    }
 
    template <typename Pusher>
    void applyFractionalStep(const Pusher& pusher, double tau, double pathLengthTarget) {
        refPartR_m /= (Physics::c * 2 * tau);
        pusher.push(refPartR_m, refPartP_m, tau);
        refPartR_m *= (Physics::c * 2 * tau);
 
        sPos_m = pathLengthTarget;
        toLabTrafo_m.transformFrom(refPartR_m);
        Vector_t<double, 3> R = refPartR_m;
        toLabTrafo_m.rotateFrom(refPartP_m);
        Vector_t<double, 3> P = refPartP_m;
        CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t<double, 3>(0, 0, 1)));
        toLabTrafo_m = toLabTrafo_m * update.inverted();
    }
 
    void scaleDtByCharge() {
        auto dtView       = dt.getView();
        const size_type n = this->getLocalNum();
        if (n == 0) {
            return;
        }
 
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto qView = QAttr.getView();
            Kokkos::parallel_for(
                    "ParticleContainer::scaleDtByCharge(attrs)", n,
                    KOKKOS_LAMBDA(const size_type i) { dtView(i) *= qView(i); });
        } else {
            auto QView = QView_m;
            Kokkos::parallel_for(
                    "ParticleContainer::scaleDtByCharge(single)", n,
                    KOKKOS_LAMBDA(const size_type i) { dtView(i) *= QView(0); });
        }
        Kokkos::fence();
    }
 
    void unscaleDtByCharge() {
        auto dtView       = dt.getView();
        const size_type n = this->getLocalNum();
        if (n == 0) {
            return;
        }
 
        if (qmStorageMode_m == QMStorageMode::Attributes) {
            auto qView = QAttr.getView();
            Kokkos::parallel_for(
                    "ParticleContainer::unscaleDtByCharge(attrs)", n,
                    KOKKOS_LAMBDA(const size_type i) { dtView(i) /= qView(i); });
        } else {
            auto QView = QView_m;
            Kokkos::parallel_for(
                    "ParticleContainer::unscaleDtByCharge(single)", n,
                    KOKKOS_LAMBDA(const size_type i) { dtView(i) /= QView(0); });
        }
        Kokkos::fence();
    }
 
    void switchToUnitlessPositions() {
        if (containerState_m->unitlessPositions) {
            throw OpalException(
                    "ParticleContainer::switchToUnitlessPositions",
                    "ParticleContainer is already in unitless positions!");
        }
        auto Rview             = this->R.getView();
        auto dtview            = this->dt.getView();
        const size_type nLocal = this->getLocalNum();
        Kokkos::parallel_for(
                "ParticleContainer::switchToUnitlessPositions", nLocal,
                KOKKOS_LAMBDA(const size_type i) { Rview(i) *= 1.0 / (Physics::c * dtview(i)); });
        Kokkos::fence();
        containerState_m->setUnitlessPositions(true);
    }
 
    void switchOffUnitlessPositions() {
        if (!containerState_m->unitlessPositions) {
            throw OpalException(
                    "ParticleContainer::switchOffUnitlessPositions",
                    "ParticleContainer is already in physical positions!");
        }
        auto Rview             = this->R.getView();
        auto dtview            = this->dt.getView();
        const size_type nLocal = this->getLocalNum();
        Kokkos::parallel_for(
                "ParticleContainer::switchOffUnitlessPositions", nLocal,
                KOKKOS_LAMBDA(const size_type i) { Rview(i) *= Physics::c * dtview(i); });
        Kokkos::fence();
        containerState_m->setUnitlessPositions(false);
    }
    QMStorageMode getQMStorageMode() const { return qmStorageMode_m; }
 
    void setGlobalProcesses(std::vector<std::unique_ptr<GlobalProcess>> processes) {
        globalProcesses_m = std::move(processes);
    }
 
    const std::vector<std::unique_ptr<GlobalProcess>>& getGlobalProcesses() const {
        return globalProcesses_m;
    }
 
    size_type markParticlesOutside(double sigmasAway) {
        size_type nLocal = this->getLocalNum();
 
        if (nLocal == 0 && this->getTotalNum() == 0) return 0;
        if (sigmasAway <= 0.0) return 0;
 
        // Force fresh moments: cached values from bunchUpdate may not reflect the exact
        // particle state at this point (emission/migration ordering can shift R between
        // the last bunchUpdate and this call). Safety-critical deletion always recomputes.
        markMomentsDirty();
        updateMoments();
 
        Vector_t<double, Dim> meanR = getMeanR();
        Vector_t<double, Dim> rmsR  = getRmsR();
 
        double lb0 = meanR[0] - sigmasAway * rmsR[0];
        double lb1 = meanR[1] - sigmasAway * rmsR[1];
        double lb2 = meanR[2] - sigmasAway * rmsR[2];
        double ub0 = meanR[0] + sigmasAway * rmsR[0];
        double ub1 = meanR[1] + sigmasAway * rmsR[1];
        double ub2 = meanR[2] + sigmasAway * rmsR[2];
 
        auto invalid = InvalidMask.getView();
        auto Rview   = this->R.getView();
 
        size_type localMarkedNum = 0;
        Kokkos::parallel_reduce(
                "ParticleContainer::markParticlesOutside", nLocal,
                KOKKOS_LAMBDA(const size_type i, size_type& count) {
                    bool outside = (Rview(i)[0] < lb0 || Rview(i)[0] > ub0)
                                   || (Rview(i)[1] < lb1 || Rview(i)[1] > ub1)
                                   || (Rview(i)[2] < lb2 || Rview(i)[2] > ub2);
                    const bool newlyMarked = outside && !invalid(i);
                    invalid(i)             = invalid(i) || outside;
                    count += newlyMarked ? 1 : 0;
                },
                localMarkedNum);
        Kokkos::fence();  // not needed, but also doesn't hurt
 
        size_type globalMarkedNum = 0;
        ippl::Comm->allreduce(localMarkedNum, globalMarkedNum, 1, std::plus<size_type>());
        return globalMarkedNum;
    }
 
    void createParticles(size_type numParticles) {
        Inform m("ParticleContainer::createParticles");
 
        // Total allocated capacity of the underlying view
        size_type oldCapacity = this->R.size();
        size_type oldLocalNum = this->getLocalNum();
        this->create(numParticles, true);  // non_destructive = true
        size_type newCapacity = this->R.size();
 
        // Initialize newly active slots.  Emitted particles are created after
        // the per-step field reset, so E/B must be explicitly cleared here.
        auto invalid = InvalidMask.getView();
        auto dtView  = dt.getView();
        auto phiView = Phi.getView();
        auto binView = Bin.getView();
        auto eView   = E.getView();
        auto bView   = B.getView();
        Kokkos::parallel_for(
                "ParticleContainer::createParticles::initializeNewSlots", numParticles,
                KOKKOS_LAMBDA(const size_type i) {
                    const size_type idx = oldLocalNum + i;
                    invalid(idx)        = false;
                    dtView(idx)         = 0.0;
                    phiView(idx)        = 0.0;
                    binView(idx)        = 0;
                    eView(idx)          = Vector_t<double, Dim>(0.0);
                    bView(idx)          = Vector_t<double, Dim>(0.0);
                });
        Kokkos::fence();
 
        // Pretty print numParticles, newCapacity and new totalNum + localNum after creation
        constexpr int labelWidth = 32;
        m << level4 << std::left << std::setw(labelWidth) << "Requested creation:" << numParticles
          << " particles" << endl;
        m << level4 << std::setw(labelWidth) << "New total number:" << this->getTotalNum()
          << " (local: " << this->getLocalNum() << ")" << endl;
        m << level4 << std::setw(labelWidth) << "Underlying view capacity:" << newCapacity << endl;
 
        if (newCapacity != oldCapacity) {
            m << level1
              << "WARNING: createParticles triggered a reallocation of the underlying particle "
                 "views! This can be a costly operation. To avoid this, consider increasing "
                 "preallocation (BEAM::NALLOC) or the overallocation factor."
              << endl;
        }
    }
 
    void allocateParticles(size_type numParticles) {
        Inform m("ParticleContainer::allocateParticles");
 
        // Total allocated capacity of the underlying view
        size_type oldCapacity = this->R.size();
        if (oldCapacity != 0) {
            throw OpalException(
                    "ParticleContainer::allocateParticles",
                    "Underlying views already allocated. This function is meant to be called on an "
                    "empty container, since it is destructive on existing particles. If you want "
                    "to create particles without deallocating existing ones, use createParticles() "
                    "instead.");
        }
 
        this->alloc(numParticles);  // alloc is always destructive
 
        m << level4 << std::left << std::setw(32) << "Requested allocation:" << numParticles
          << " particles" << endl;
        m << level4 << std::setw(32) << "Size of underlying view:" << this->R.size() << endl;
    }
 
    size_type deleteInvalidParticles() {
        Inform m("ParticleContainer::deleteInvalidParticles");
 
        const size_type nLocal = this->getLocalNum();
        auto invalid           = InvalidMask.getView();
        if (invalid.extent(0) < nLocal) {
            throw OpalException(
                    "ParticleContainer::deleteInvalidParticles",
                    "InvalidMask extent (" + std::to_string(invalid.extent(0))
                            + ") is smaller than local particle count (" + std::to_string(nLocal)
                            + ").");
        }
 
        size_type localDestroyNum = 0;
        Kokkos::parallel_reduce(
                "ParticleContainer::deleteInvalidParticles::count", nLocal,
                KOKKOS_LAMBDA(const size_type i, size_type& count) { count += invalid(i) ? 1 : 0; },
                localDestroyNum);
        Kokkos::fence();
 
        size_type globalDestroyNum = 0;
        ippl::Comm->allreduce(localDestroyNum, globalDestroyNum, 1, std::plus<size_type>());
 
        // This is a collective call! All ranks must execute it. Note that this is safe with the
        // current IPPL implementation: ParticleBase::internalDestroy() first consumes the invalid
        // view to build internal deleteIndex_m and keepIndex_m arrays. Meaning, the fact that
        // destroy edits the InvalidMask while using it is not a problem.
        Base::destroy(invalid, localDestroyNum);
 
        // Mark moments dirty only if there were changes globally.
        if (globalDestroyNum > 0) {
            markMomentsDirty();
        }
 
        // Reset particle mask after deletion.
        InvalidMask = false;
        Kokkos::fence();
 
        constexpr int labelWidth = 32;
        m << level4 << std::left << std::setw(labelWidth)
          << "Requested destruction:" << localDestroyNum << " particles" << endl
          << std::setw(labelWidth) << "New total number:" << this->getTotalNum()
          << " (local: " << this->getLocalNum() << ")" << endl;
 
        return globalDestroyNum;
    }
 
private:
    void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); }
 
    PLayout_t<T, Dim> pl_m;
 
    QMStorageMode qmStorageMode_m = QMStorageMode::SingleValue;
 
    DistributionMoments distMoments_m;
 
    std::shared_ptr<BunchStateHandler::ContainerState> containerState_m;
 
    // Single shared scalar mode stored as a length-1 Kokkos view.
    qm_view_type QView_m;
    qm_view_type MView_m;
 
    // Per-particle attributes mode
    ippl::ParticleAttrib<double> QAttr;
    ippl::ParticleAttrib<double> MAttr;
 
    // Whether the spin attribute is registered on this container.
    bool spinEnabled_m = false;
 
    // Reference particle information
    Vector_t<double, Dim> refPartR_m;
    Vector_t<double, Dim> refPartP_m;
 
    // Global to local quaternion
    Quaternion_t globalToLocalQuaternion_m;
 
    // Reference particle to lab transformation
    CoordinateSystemTrafo toLabTrafo_m;
 
    // Particle reference data (!= reference particle)
    const PartData* reference_m = nullptr;
 
    // Distance along the beamline
    double sPos_m = 0.0;
 
    std::vector<std::unique_ptr<GlobalProcess>> globalProcesses_m;
};
 
#endif