DistributionMoments.h

Go to the documentation of this file.

//
// Class DistributionMoments
//   Computes the statistics of particle distributions.
//
// Copyright (c) 2025, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef DISTRIBUTIONMOMENTS_H
#define DISTRIBUTIONMOMENTS_H
 
#include "Algorithms/Matrix.h"
#include "Ippl.h"
#include "PartBunch/BunchStateHandler.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
 
#include <Kokkos_Core.hpp>
#include <vector>
 
template <typename T, unsigned Dim = 3>
using Vector_t = ippl::Vector<T, Dim>;
 
typedef typename std::pair<Vector_t<double, 3>, Vector_t<double, 3>> VectorPair_t;
 
#include "Algorithms/Matrix.h"
 
class OpalParticle;
 
class DistributionMoments {
public:
    DistributionMoments();
 
    void setContainerState(std::weak_ptr<BunchStateHandler::ContainerState> containerState) {
        containerState_m = std::move(containerState);
    }
 
    void setEnergyReferenceMass(double referenceMassGeV, bool rescaleToReference = true) {
        referenceMassGeV_m   = referenceMassGeV;
        rescaleToReference_m = rescaleToReference;
    }
    bool getRescaleEnergyToReference() const { return rescaleToReference_m; }
    double getEnergyReferenceMassGeV() const { return referenceMassGeV_m; }
    void compute(
            const std::vector<OpalParticle>::const_iterator&,
            const std::vector<OpalParticle>::const_iterator&);
    void computeMoments(
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Rview,
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Pview,
            ippl::ParticleAttrib<double>::view_type Mview, size_t Np, size_t Nlocal);
    void computeMinMaxPosition(
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Rview, size_t Nlcoal);
    void computeMeanKineticEnergy();
    void computeDebyeLength(
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Pview, size_t Np, size_t Nlocal,
            double density);
    void computePlasmaParameter(double);
 
    Vector_t<double, 3> getMeanPosition() const;
    Vector_t<double, 3> getStandardDeviationPosition() const;
    Vector_t<double, 3> getMeanMomentum() const;
    Vector_t<double, 3> getStandardDeviationMomentum() const;
    Vector_t<double, 3> getNormalizedEmittance() const;
    Vector_t<double, 3> getGeometricEmittance() const;
    Vector_t<double, 3> getStandardDeviationRP() const;
    Vector_t<double, 3> getHalo() const;
    Vector_t<double, 3> getMinPosition() const;
    Vector_t<double, 3> getMaxPosition() const;
    Vector_t<double, 3> getMaxR() const;
 
    Vector_t<double, 3> get68Percentile() const;
    Vector_t<double, 3> get95Percentile() const;
    Vector_t<double, 3> get99Percentile() const;
    Vector_t<double, 3> get99_99Percentile() const;
    Vector_t<double, 3> getNormalizedEmittance68Percentile() const;
    Vector_t<double, 3> getNormalizedEmittance95Percentile() const;
    Vector_t<double, 3> getNormalizedEmittance99Percentile() const;
    Vector_t<double, 3> getNormalizedEmittance99_99Percentile() const;
 
    double getMeanTime() const;
    double getStdTime() const;
    double getMeanGamma() const;
    double getMeanGammaZ() const;
    double getMeanKineticEnergy() const;
    double getTemperature() const;
    double getDebyeLength() const;
    double getPlasmaParameter() const;
    double getStdKineticEnergy() const;
    double getDx() const;
    double getDDx() const;
    double getDy() const;
    double getDDy() const;
    Vector_t<double, 6> getMeans() const;
    Vector_t<double, 6> getCentroid() const;
    matrix6x6_t getMoments6x6() const;
    double getTotalCharge() const;
    double getTotalMass() const;
    double getTotalNumParticles() const;
 
    void computeMeans(
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Rview,
            ippl::ParticleAttrib<Vector_t<double, 3>>::view_type Pview,
            ippl::ParticleAttrib<double>::view_type Mview, size_t Np, size_t Nlocal);
 
private:
    bool isParticleExcluded(const OpalParticle&) const;
 
    // template <class InputIt>
    // void computeMeans(const InputIt&, const InputIt&);
    //  template <class InputIt>
    //  void computeStatistics(const InputIt&, const InputIt&);
    template <class InputIt>
    void computePercentiles(const InputIt&, const InputIt&);
    using iterator_t = std::vector<Vector_t<double, 2>>::const_iterator;
    std::pair<double, iterator_t> determinePercentilesDetail(
            const iterator_t& begin, const iterator_t& end,
            const std::vector<int>& globalAccumulatedHistogram,
            const std::vector<int>& localAccumulatedHistogram, unsigned int dimension,
            int numRequiredParticles) const;
    double computeNormalizedEmittance(const iterator_t& begin, const iterator_t& end) const;
 
    void fillMembers(std::vector<double>&);
 
    void reset();
 
    void resetPlasmaParameters();
 
    std::weak_ptr<BunchStateHandler::ContainerState> containerState_m;
 
    Vector_t<double, 3> meanR_m;
    Vector_t<double, 3> meanP_m;
    Vector_t<double, 3> stdR_m;
    Vector_t<double, 3> stdP_m;
    Vector_t<double, 3> stdRP_m;
    Vector_t<double, 3> normalizedEps_m;
    Vector_t<double, 3> geometricEps_m;
    Vector_t<double, 3> halo_m;
    Vector_t<double, 3> maxR_m;
    Vector_t<double, 3> minR_m;
    Vector_t<double, 3> sixtyEightPercentile_m;
    Vector_t<double, 3> normalizedEps68Percentile_m;
    Vector_t<double, 3> ninetyFivePercentile_m;
    Vector_t<double, 3> normalizedEps95Percentile_m;
    Vector_t<double, 3> ninetyNinePercentile_m;
    Vector_t<double, 3> normalizedEps99Percentile_m;
    Vector_t<double, 3> ninetyNine_NinetyNinePercentile_m;
    Vector_t<double, 3> normalizedEps99_99Percentile_m;
 
    double meanTime_m;
    double stdTime_m;
    double meanKineticEnergy_m;
    double temperature_m;
    double debyeLength_m;
    double plasmaParameter_m;
    double stdKineticEnergy_m;
    double meanGamma_m;
    double meanGammaZ_m;
 
    // If enabled, compute kinetic energy moments with referenceMassGeV_m instead of Mview(k).
    bool rescaleToReference_m = false;
    double referenceMassGeV_m = 0.0;
 
    Vector_t<double, 6> centroid_m;
    Vector_t<double, 6> means_m;
    matrix6x6_t moments_m;
    matrix6x6_t notCentMoments_m;
 
    double totalCharge_m;
    double totalMass_m;
    unsigned int totalNumParticles_m;
 
    static const double percentileOneSigmaNormalDist_m;
    static const double percentileTwoSigmasNormalDist_m;
    static const double percentileThreeSigmasNormalDist_m;
    static const double percentileFourSigmasNormalDist_m;
};
 
inline Vector_t<double, 3> DistributionMoments::getMeanPosition() const { return meanR_m; }
 
inline Vector_t<double, 3> DistributionMoments::getStandardDeviationPosition() const {
    return stdR_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getMeanMomentum() const { return meanP_m; }
 
inline Vector_t<double, 3> DistributionMoments::getStandardDeviationMomentum() const {
    return stdP_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getNormalizedEmittance() const {
    return normalizedEps_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getGeometricEmittance() const {
    return geometricEps_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getStandardDeviationRP() const { return stdRP_m; }
 
inline Vector_t<double, 3> DistributionMoments::getHalo() const { return halo_m; }
 
inline Vector_t<double, 3> DistributionMoments::getMinPosition() const { return minR_m; }
 
inline Vector_t<double, 3> DistributionMoments::getMaxPosition() const { return maxR_m; }
 
inline double DistributionMoments::getMeanTime() const { return meanTime_m; }
 
inline double DistributionMoments::getStdTime() const { return stdTime_m; }
 
inline double DistributionMoments::getMeanGamma() const { return meanGamma_m; }
 
inline double DistributionMoments::getMeanGammaZ() const { return meanGammaZ_m; }
 
inline double DistributionMoments::getMeanKineticEnergy() const { return meanKineticEnergy_m; }
 
// Compute and return the value of temperature in K
inline double DistributionMoments::getTemperature() const {
    return (temperature_m / (Physics::kB * Units::eV2kg * Physics::c * Physics::c));
}
inline double DistributionMoments::getDebyeLength() const { return debyeLength_m; }
inline double DistributionMoments::getPlasmaParameter() const { return plasmaParameter_m; }
 
inline double DistributionMoments::getStdKineticEnergy() const { return stdKineticEnergy_m; }
 
inline double DistributionMoments::getDx() const { return moments_m(0, 5); }
 
inline double DistributionMoments::getDDx() const { return moments_m(1, 5); }
 
inline double DistributionMoments::getDy() const { return moments_m(2, 5); }
 
inline double DistributionMoments::getDDy() const { return moments_m(3, 5); }
 
inline Vector_t<double, 6> DistributionMoments::getCentroid() const { return centroid_m; }
 
inline Vector_t<double, 6> DistributionMoments::getMeans() const { return means_m; }
 
inline matrix6x6_t DistributionMoments::getMoments6x6() const { return moments_m; }
 
inline double DistributionMoments::getTotalCharge() const { return totalCharge_m; }
 
inline double DistributionMoments::getTotalMass() const { return totalMass_m; }
 
inline double DistributionMoments::getTotalNumParticles() const { return totalNumParticles_m; }
 
inline Vector_t<double, 3> DistributionMoments::get68Percentile() const {
    return sixtyEightPercentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getNormalizedEmittance68Percentile() const {
    return normalizedEps68Percentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::get95Percentile() const {
    return ninetyFivePercentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getNormalizedEmittance95Percentile() const {
    return normalizedEps95Percentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::get99Percentile() const {
    return ninetyNinePercentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getNormalizedEmittance99Percentile() const {
    return normalizedEps99Percentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::get99_99Percentile() const {
    return ninetyNine_NinetyNinePercentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getNormalizedEmittance99_99Percentile() const {
    return normalizedEps99_99Percentile_m;
}
 
inline Vector_t<double, 3> DistributionMoments::getMaxR() const {
    Vector_t<double, 3> maxDistance;
    for (unsigned int i = 0; i < 3; ++i) {
        maxDistance[i] = std::max(std::abs(maxR_m[i]), std::abs(minR_m[i]));
    }
    return maxDistance;
}
 
// ---------------------------------------------------------------------------
// Custom Kokkos-reducible types for parallel_reduce.
// All used with Kokkos::Sum (which calls operator+= and reduction_identity::sum()).
// See
// https://kokkos.org/kokkos-core-wiki/ProgrammingGuide/Custom-Reductions-Built-In-Reducers-with-Custom-Scalar-Types.html
// ---------------------------------------------------------------------------
 
// Element-wise sum of N doubles.
template <int N>
struct SumArray {
    double data[N];
 
    KOKKOS_INLINE_FUNCTION
    SumArray() {
        for (int i = 0; i < N; ++i)
            data[i] = 0.0;
    }
 
    KOKKOS_INLINE_FUNCTION
    SumArray& operator+=(const SumArray& src) {
        for (int i = 0; i < N; ++i)
            data[i] += src.data[i];
        return *this;
    }
};
 
// Element-wise maximum of N doubles.
template <int N>
struct MaxArray {
    double data[N];
 
    KOKKOS_INLINE_FUNCTION
    MaxArray() {
        for (int i = 0; i < N; ++i)
            data[i] = Kokkos::reduction_identity<double>::max();
    }
 
    KOKKOS_INLINE_FUNCTION
    MaxArray& operator+=(const MaxArray& src) {
        for (int i = 0; i < N; ++i)
            if (src.data[i] > data[i]) data[i] = src.data[i];
        return *this;
    }
};
 
// Element-wise minimum of N doubles.
template <int N>
struct MinArray {
    double data[N];
 
    KOKKOS_INLINE_FUNCTION
    MinArray() {
        for (int i = 0; i < N; ++i)
            data[i] = Kokkos::reduction_identity<double>::min();
    }
 
    KOKKOS_INLINE_FUNCTION
    MinArray& operator+=(const MinArray& src) {
        for (int i = 0; i < N; ++i)
            if (src.data[i] < data[i]) data[i] = src.data[i];
        return *this;
    }
};
 
// Element-wise sum of a 6x6 matrix of doubles.
struct SumMatrix6x6 {
    double data[6][6];
 
    KOKKOS_INLINE_FUNCTION
    SumMatrix6x6() {
        for (int i = 0; i < 6; ++i)
            for (int j = 0; j < 6; ++j)
                data[i][j] = 0.0;
    }
 
    KOKKOS_INLINE_FUNCTION
    SumMatrix6x6& operator+=(const SumMatrix6x6& src) {
        for (int i = 0; i < 6; ++i)
            for (int j = 0; j < 6; ++j)
                data[i][j] += src.data[i][j];
        return *this;
    }
};
 
namespace Kokkos {
    template <int N>
    struct reduction_identity<SumArray<N>> {
        KOKKOS_FORCEINLINE_FUNCTION static SumArray<N> sum() { return SumArray<N>(); }
    };
    template <int N>
    struct reduction_identity<MaxArray<N>> {
        KOKKOS_FORCEINLINE_FUNCTION static MaxArray<N> sum() { return MaxArray<N>(); }
    };
    template <int N>
    struct reduction_identity<MinArray<N>> {
        KOKKOS_FORCEINLINE_FUNCTION static MinArray<N> sum() { return MinArray<N>(); }
    };
    template <>
    struct reduction_identity<SumMatrix6x6> {
        KOKKOS_FORCEINLINE_FUNCTION static SumMatrix6x6 sum() { return SumMatrix6x6(); }
    };
}  // namespace Kokkos
 
// ---------------------------------------------------------------------------
 
#endif