29#include <gsl/gsl_histogram.h>
31#include <boost/numeric/conversion/cast.hpp>
52 const std::vector<OpalParticle>::const_iterator &last)
57template<
class InputIt>
60 unsigned int localNum = last - first;
61 std::vector<double> localStatistics(10);
62 std::vector<double> maxima(6, std::numeric_limits<double>::lowest());
63 for (InputIt it = first; it != last; ++ it) {
69 for (
unsigned int i = 0; i < 3; ++ i) {
70 maxima[2 * i] =
std::max(maxima[2 * i], particle[2 * i]);
71 maxima[2 * i + 1] =
std::max(maxima[2 * i + 1], -particle[2 * i]);
75 for (
unsigned int i = 0; i < 6; ++ i, ++ l) {
76 localStatistics[l] += particle[i];
80 localStatistics[l++] += particle.
getTime();
83 double eKin = (gamma - 1.0) * particle.
getMass();
84 localStatistics[l++] += eKin;
85 localStatistics[l++] += gamma;
87 localStatistics.back() = localNum;
89 allreduce(localStatistics.data(), localStatistics.size(), std::plus<double>());
90 allreduce(maxima.data(), 6, std::greater<double>());
99 for (
unsigned int i = 0; i < 3; ++ i) {
102 maxR_m(i) = maxima[2 * i];
103 minR_m(i) = -maxima[2 * i + 1];
106 meanTime_m = localStatistics[l++] * perParticle;
122template<
class InputIt>
130 std::vector<double> localStatistics(37);
131 for (InputIt it = first; it != last; ++ it) {
137 for (
unsigned int i = 0; i < 6; ++ i) {
139 for (
unsigned int j = 0; j <= i; ++ j, ++ l) {
140 localStatistics[l] += particle[i] * particle[j];
146 for (
unsigned int i = 0; i < 3; ++ i, l += 2) {
147 double r2 =
std::pow(particle[i], 2);
148 localStatistics[l] += r2 * particle[i];
149 localStatistics[l + 1] += r2 * r2;
154 localStatistics[l++] += particle.
getCharge();
155 localStatistics[l++] += particle.
getMass();
158 allreduce(localStatistics.data(), localStatistics.size(), std::plus<double>());
166template<
class InputIt>
172 std::vector<gsl_histogram*> histograms(3);
179 for (
unsigned int d = 0; d < 3; ++ d) {
181 histograms[d] = gsl_histogram_alloc(numBins);
182 gsl_histogram_set_ranges_uniform(histograms[d], 0.0, maxR(d));
184 for (InputIt it = first; it != last; ++ it) {
186 for (
unsigned int d = 0; d < 3; ++ d) {
187 gsl_histogram_increment(histograms[d],
std::abs(particle[2 * d] -
meanR_m[d]));
191 std::vector<int> localHistogramValues(3 * (numBins + 1)), globalHistogramValues(3 * (numBins + 1));
192 for (
unsigned int d = 0; d < 3; ++ d) {
194 size_t accumulated = 0;
195 std::generate(localHistogramValues.begin() + d * (numBins + 1) + 1,
196 localHistogramValues.begin() + (d + 1) * (numBins + 1),
197 [&histograms,&d,&j,&accumulated](){
198 accumulated += gsl_histogram_get(histograms[d], j ++);
199 return accumulated;});
201 gsl_histogram_free(histograms[d]);
204 allreduce(localHistogramValues.data(), globalHistogramValues.data(), 3 * (numBins + 1), std::plus<int>());
211 for (
int d = 0; d < 3; ++ d) {
212 unsigned int localNum = last - first, current = 0;
213 std::vector<Vektor<double, 2>> oneDPhaseSpace(localNum);
214 for (InputIt it = first; it != last; ++ it, ++ current) {
216 oneDPhaseSpace[current](0) = particle[2 * d];
217 oneDPhaseSpace[current](1) = particle[2 * d + 1];
219 std::sort(oneDPhaseSpace.begin(), oneDPhaseSpace.end(),
221 return std::abs(left[0] - meanR_m[d]) < std::abs(right[0] - meanR_m[d]);
224 iterator_t endSixtyEight, endNinetyFive, endNinetyNine, endNinetyNine_NinetyNine;
225 endSixtyEight = endNinetyFive = endNinetyNine = endNinetyNine_NinetyNine = oneDPhaseSpace.end();
229 oneDPhaseSpace.end(),
230 globalHistogramValues,
231 localHistogramValues,
236 oneDPhaseSpace.end(),
237 globalHistogramValues,
238 localHistogramValues,
243 oneDPhaseSpace.end(),
244 globalHistogramValues,
245 localHistogramValues,
250 oneDPhaseSpace.end(),
251 globalHistogramValues,
252 localHistogramValues,
253 d, numParticles99_99);
291std::pair<double, DistributionMoments::iterator_t>
294 const std::vector<int>& globalAccumulatedHistogram,
295 const std::vector<int>& localAccumulatedHistogram,
296 unsigned int dimension,
297 int numRequiredParticles)
const {
298 unsigned int numBins = globalAccumulatedHistogram.size() / 3;
299 double percentile = 0.0;
301 for (
unsigned int i = 1; i < numBins; ++ i) {
302 unsigned int idx = dimension * numBins + i;
303 if (globalAccumulatedHistogram[idx] > numRequiredParticles) {
306 unsigned int numMissingParticles = numRequiredParticles - globalAccumulatedHistogram[idx - 1];
307 unsigned int shift = 2;
308 while (numMissingParticles == 0) {
309 beginBin =
begin + localAccumulatedHistogram[idx - shift];
310 numMissingParticles = numRequiredParticles - globalAccumulatedHistogram[idx - shift];
314 std::vector<unsigned int> numParticlesInBin(
Ippl::getNodes() + 1);
315 numParticlesInBin[
Ippl::myNode() + 1] = endBin - beginBin;
317 std::partial_sum(numParticlesInBin.begin(), numParticlesInBin.end(), numParticlesInBin.begin());
319 std::vector<double> positions(numParticlesInBin.back());
320 std::transform(beginBin, endBin, positions.begin() + numParticlesInBin[
Ippl::myNode()],
322 { return std::abs(particle[0] - meanR_m[dimension]); });
323 allreduce(&(positions[0]), positions.size(), std::plus<double>());
324 std::sort(positions.begin(), positions.end());
326 percentile = (*(positions.begin() + numMissingParticles - 1)
327 + *(positions.begin() + numMissingParticles)) / 2;
328 for (
iterator_t it = beginBin; it != endBin; ++ it) {
330 return std::make_pair(percentile, it);
333 endPercentile = endBin;
337 return std::make_pair(percentile, endPercentile);
342 double localStatistics[] = {0.0, 0.0, 0.0, 0.0};
346 localStatistics[1] += rp(0);
347 localStatistics[2] += rp(1);
348 localStatistics[3] += rp(0) * rp(1);
350 allreduce(&(localStatistics[0]), 4, std::plus<double>());
352 double numParticles = localStatistics[0];
353 double perParticle = 1 / localStatistics[0];
354 double meanR = localStatistics[1] * perParticle;
355 double meanP = localStatistics[2] * perParticle;
356 double RP = localStatistics[3] * perParticle;
358 localStatistics[0] = 0.0;
359 localStatistics[1] = 0.0;
362 localStatistics[0] +=
std::pow(rp(0) - meanR, 2);
363 localStatistics[1] +=
std::pow(rp(1) - meanP, 2);
365 allreduce(&(localStatistics[0]), 2, std::plus<double>());
367 double stdR =
std::sqrt(localStatistics[0] / numParticles);
368 double stdP =
std::sqrt(localStatistics[1] / numParticles);
369 double sumRP = RP - meanR * meanP;
373 return normalizedEps;
381 for (; l < 6; l += 2) {
386 for (
unsigned int i = 0; i < 6; ++ i) {
387 for (
unsigned int j = 0; j <= i; ++ j, ++ l) {
388 moments_m(i, j) = localMoments[l] * perParticle;
393 stdTime_m = localMoments[l++] * perParticle;
395 for (
unsigned int i = 0; i < 3; ++ i, l += 2) {
398 double w3 = localMoments[l] * perParticle;
399 double w4 = localMoments[l + 1] * perParticle;
402 halo_m(i) = (w4 + w1 * (-4 * w3 + 3 * w1 * (tmp + w2))) / tmp;
411 for (
unsigned int i = 0; i < 3; ++ i) {
424 double data[] = {0.0, 0.0};
438 double avgVel[3]={0.0,0.0,0.0};
442 for(
unsigned i = 0; i < 3; i++) {
443 avgVel[i] += ((particle_r.getP()[i] *
Physics::c)/
448 allreduce(avgVel, 3, std::plus<double>());
450 const double N =
static_cast<double>(bunch_r.
getTotalNum());
451 for(
unsigned i = 0; i < 3; i++) {
452 avgVel[i]= avgVel[i]/N;
455 double tempAvg = 0.0;
458 for(
unsigned i = 0; i < 3; i++) {
463 allreduce(tempAvg, 1, std::plus<double>());
526 && particle.
getId() == 0;
FTps< T, N > erf(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Error function.
PartBunchBase< T, Dim >::ConstIterator end(PartBunchBase< T, Dim > const &bunch)
PartBunchBase< T, Dim >::ConstIterator begin(PartBunchBase< T, Dim > const &bunch)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Tps< T > sqrt(const Tps< T > &x)
Square root.
void allreduce(const T *input, T *output, int count, Op op)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
PETE_TUTree< FnLog10, typename T::PETE_Expr_t > log10(const PETE_Expr< T > &l)
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
constexpr double epsilon_0
The permittivity of vacuum in As/Vm.
constexpr double q_e
The elementary charge in As.
constexpr double m_e
The electron rest mass in GeV.
constexpr double c
The velocity of light in m/s.
constexpr double pi
The value of.
double haloShift
The constant parameter C to shift halo, by < w^4 > / < w^2 > ^2 - C (w=x,y,z)
bool amr
Enable AMR if true.
double getKineticEnergy(Vector_t p, double mass)
double getGamma(Vector_t p)
size_t getLocalNum() const
size_t getTotalNum() const
ConstIterator end() const
ConstIterator begin() const
static OpalData * getInstance()
double meanKineticEnergy_m
Vector_t ninetyNinePercentile_m
void computeDebyeLength(PartBunchBase< double, 3 > const &, double)
void computePlasmaParameter(double)
FMatrix< double, 6, 6 > moments_m
static const double percentileFourSigmasNormalDist_m
void computeMeans(const InputIt &, const InputIt &)
Vector_t normalizedEps99_99Percentile_m
void computeMeanKineticEnergy(PartBunchBase< double, 3 > const &)
unsigned int totalNumParticles_m
void compute(const std::vector< OpalParticle >::const_iterator &, const std::vector< OpalParticle >::const_iterator &)
Vector_t normalizedEps68Percentile_m
Vector_t normalizedEps95Percentile_m
std::vector< Vektor< double, 2 > >::const_iterator iterator_t
Vector_t ninetyFivePercentile_m
void computeStatistics(const InputIt &, const InputIt &)
double computeNormalizedEmittance(const iterator_t &begin, const iterator_t &end) const
Vector_t normalizedEps99Percentile_m
static const double percentileThreeSigmasNormalDist_m
void computePercentiles(const InputIt &, const InputIt &)
static const double percentileTwoSigmasNormalDist_m
bool isParticleExcluded(const OpalParticle &) const
Vector_t sixtyEightPercentile_m
Vector_t ninetyNine_NinetyNinePercentile_m
void fillMembers(std::vector< double > const &)
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
static const double percentileOneSigmaNormalDist_m
void resetPlasmaParameters()
double stdKineticEnergy_m
const Vector_t & getP() const
Get momentum.
double getCharge() const
Get charge in Coulomb.
double getMass() const
Get mass in GeV/c^2.
double getTime() const
Get time.
int64_t getId() const
Get the id of the particle.
Vektor< double, 3 > Vector_t