Commit 8631137e authored by kraus's avatar kraus
Browse files

add initial version

parent 131ad823
......@@ -3,6 +3,7 @@ set (_SRCS
AbstractTracker.cpp
CoordinateSystemTrafo.cpp
DefaultVisitor.cpp
DistributionMoments.cpp
Flagger.cpp
PartBunch.cpp
PartBins.cpp
......@@ -30,6 +31,7 @@ set (HDRS
AbstractTracker.h
CoordinateSystemTrafo.h
DefaultVisitor.h
DistributionMoments.h
Flagger.h
ListElem.h
PartBinsCyc.h
......
#include "DistributionMoments.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "Message/GlobalComm.h"
DistributionMoments::DistributionMoments():
meanR_m(0.0),
meanP_m(0.0),
stdR_m(0.0),
stdP_m(0.0),
normalizedEps_m(0.0),
geometricEps_m(0.0),
meanKineticEnergy_m(0.0),
moments_m(),
totalNumParticles_m(0)
{
std::fill(std::begin(centroid_m), std::end(centroid_m), 0.0);
}
void DistributionMoments::compute(std::vector<OpalParticle> const& particles)
{
reset();
unsigned int localNumParticles = particles.size();
allreduce(&localNumParticles, &totalNumParticles_m, 1, std::plus<unsigned int>());
if (totalNumParticles_m == 0) {
return;
}
computeMoments(particles);
Vector_t squaredEps, fac, squaredSumR, squaredSumP, sumRP;
double perParticle = 1.0 / totalNumParticles_m;
for (unsigned int i = 0; i < 3; ++ i) {
meanR_m(i) = centroid_m[2 * i] * perParticle;
meanP_m(i) = centroid_m[2 * i + 1] * perParticle;
squaredSumR(i) = moments_m(2 * i, 2 * i) - totalNumParticles_m * std::pow(meanR_m(i), 2);
squaredSumP(i) = std::max(0.0,
moments_m(2 * i + 1, 2 * i + 1) - totalNumParticles_m * std::pow(meanP_m(i), 2));
sumRP(i) = moments_m(2 * i, 2 * i + 1) - totalNumParticles_m * meanR_m(i) * meanP_m(i);
}
squaredEps = (squaredSumR * squaredSumP - sumRP * sumRP) * std::pow(perParticle, 2);
// sumRP *= perParticle;
for (unsigned int i = 0; i < 3; ++ i) {
stdR_m(i) = std::sqrt(squaredSumR(i) * perParticle);
stdP_m(i) = std::sqrt(squaredSumP(i) * perParticle);
normalizedEps_m(i) = std::sqrt(std::max(squaredEps(i), 0.0));
double tmp = stdR_m(i) * stdP_m(i);
fac(i) = (std::abs(tmp) < 1e-10) ? 0.0: 1.0 / tmp;
}
// stdRP_m = sumRP * fac;
double betaGamma = std::sqrt(std::pow(meanGamma_m, 2) - 1.0);
geometricEps_m = normalizedEps_m / Vector_t(betaGamma);
}
void DistributionMoments::computeMoments(std::vector<OpalParticle> const& particles)
{
std::vector<double> localMoments(35);
for (OpalParticle const& particle: particles) {
unsigned int l = 6;
for (unsigned int i = 0; i < 6; ++ i) {
localMoments[i] += particle[i];
for (unsigned int j = 0; j <= i; ++ j, ++ l) {
localMoments[l] += particle[i] * particle[j];
}
}
for (unsigned int i = 0; i < 3; ++ i, ++ l) {
double r2 = std::pow(particle[i], 2);
localMoments[l] += r2 * particle[i];
localMoments[l + 3] += r2 * r2;
}
double gamma = Util::getGamma(particle.P());
localMoments[33] = (gamma - 1.0) * particle.mass();
localMoments[34] = gamma;
}
allreduce(localMoments.data(), localMoments.size(), std::plus<double>());
for (unsigned int i = 0; i < 6; ++ i) {
centroid_m[i] = localMoments[i];
}
unsigned int l = 6;
for (unsigned int i = 0; i < 6; ++ i) {
for (unsigned int j = 0; j <= i; ++ j, ++ l) {
moments_m(i, j) = localMoments[l];
moments_m(j, i) = moments_m(i, j);
}
}
double perParticle = 1.0 / totalNumParticles_m;
for (unsigned int i = 0; i < 3; ++ i) {
double w1 = centroid_m[2 * i] * perParticle;
double w2 = moments_m(2 * i , 2 * i) * perParticle;
double w3 = localMoments[i + l] * perParticle;
double w4 = localMoments[i + l + 3] * perParticle;
double tmp = w2 - std::pow(w1, 2);
halo_m(i) = (w4 + w1 * (-4 * w3 + 3 * w1 * (tmp + w2))) / tmp;
halo_m(i) -= Options::haloShift;
}
meanKineticEnergy_m = localMoments[33] * perParticle;
meanGamma_m = localMoments[34] * perParticle;
}
void DistributionMoments::reset()
{
meanR_m = 0.0;
meanP_m = 0.0;
stdR_m = 0.0;
stdP_m = 0.0;
normalizedEps_m = 0.0;
geometricEps_m = 0.0;
halo_m = 0.0;
meanKineticEnergy_m = 0.0;
std::fill(std::begin(centroid_m), std::end(centroid_m), 0.0);
moments_m = FMatrix<double, 6, 6>(0.0);
totalNumParticles_m = 0;
}
\ No newline at end of file
#ifndef DISTRIBUTIONMOMENTS_H
#define DISTRIBUTIONMOMENTS_H
#include "OpalParticle.h"
#include "FixedAlgebra/FMatrix.h"
#include <vector>
class DistributionMoments {
public:
DistributionMoments();
void compute(std::vector<OpalParticle> const&);
Vector_t getMeanPosition() const;
Vector_t getStandardDeviationPosition() const;
Vector_t getMeanMomentum() const;
Vector_t getStandardDeviationMomentum() const;
Vector_t getNormalizedEmittance() const;
Vector_t getGeometricEmittance() const;
double getMeanKineticEnergy() const;
private:
void computeMoments(std::vector<OpalParticle> const&);
void reset();
Vector_t meanR_m;
Vector_t meanP_m;
Vector_t stdR_m;
Vector_t stdP_m;
Vector_t normalizedEps_m;
Vector_t geometricEps_m;
Vector_t halo_m;
double meanKineticEnergy_m;
double meanGamma_m;
double centroid_m[6];
FMatrix<double, 6, 6> moments_m;
unsigned int totalNumParticles_m;
};
inline
Vector_t DistributionMoments::getMeanPosition() const
{
return meanR_m;
}
inline
Vector_t DistributionMoments::getStandardDeviationPosition() const
{
return stdR_m;
}
inline
Vector_t DistributionMoments::getMeanMomentum() const
{
return meanP_m;
}
inline
Vector_t DistributionMoments::getStandardDeviationMomentum() const
{
return stdP_m;
}
inline
Vector_t DistributionMoments::getNormalizedEmittance() const
{
return normalizedEps_m;
}
inline
Vector_t DistributionMoments::getGeometricEmittance() const
{
return geometricEps_m;
}
inline
double DistributionMoments::getMeanKineticEnergy() const
{
return meanKineticEnergy_m;
}
#endif
\ No newline at end of file
// ------------------------------------------------------------------------
// $RCSfile: OpalParticle.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: OpalParticle
// A OpalParticle represents the phase space coordinates of a particle.
// It can be propagated through a beamline.
// Class OpalParticle
// This class represents the canonical coordinates of a particle.
//
// ------------------------------------------------------------------------
// Class category: Algorithms
// ------------------------------------------------------------------------
// Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// $Date: 2000/03/27 09:32:33 $
// $Author: fci $
// 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/>.
//
// ------------------------------------------------------------------------
#include "Algorithms/OpalParticle.h"
// Class OpalParticle
// ------------------------------------------------------------------------
OpalParticle::OpalParticle()
{}
OpalParticle::OpalParticle
(double x, double px, double y, double py, double t, double pt)
{
phase[X] = x;
phase[PX] = px;
phase[Y] = y;
phase[PY] = py;
phase[T] = t;
phase[PT] = pt;
}
OpalParticle::OpalParticle(double x, double px,
double y, double py,
double t, double pt,
double q, double m):
R_m(x, y, t),
P_m(px, py, pt),
charge_m(q),
mass_m(m)
{}
\ No newline at end of file
#ifndef CLASSIC_OpalParticle_HH
#define CLASSIC_OpalParticle_HH
// ------------------------------------------------------------------------
// $RCSfile: OpalParticle.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: OpalParticle
// Class OpalParticle
// This class represents the canonical coordinates of a particle.
//
// ------------------------------------------------------------------------
// Class category: Algorithms
// ------------------------------------------------------------------------
// Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// $Date: 2000/03/27 09:32:33 $
// $Author: fci $
// 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 CLASSIC_OpalParticle_HH
#define CLASSIC_OpalParticle_HH
// Class OpalParticle
// ------------------------------------------------------------------------
/// OpalParticle position.
// This class represents the canonical coordinates of a particle.
// {P}
// NOTE. The order of phase space coordinates is,
// as opposed to the original CLASSIC definition:
// {CENTER}
// X = 0, PX = 1, Y = 2, PY = 3, T = 4, PT = 5.
// {/CENTER}
// The copy constructor, destructor, and assignment operator generated
// by the compiler perform the correct operation. For speed reasons
// they are not implemented.
#include "Vektor.h"
class OpalParticle {
......@@ -44,7 +29,10 @@ public:
/// Constructor.
// Construct particle with the given coordinates.
OpalParticle(double x, double px, double y, double py, double t, double pt);
OpalParticle(double x, double px,
double y, double py,
double t, double pt,
double q, double m);
OpalParticle();
......@@ -71,6 +59,12 @@ public:
/// Get reference to relative momentum error (no dimension).
double &pt();
/// Get reference to position in m
Vector_t &R();
/// Get reference to momentum
Vector_t &P();
/// Get coordinate.
// Access coordinate by index for constant particle.
double operator[](int) const;
......@@ -93,61 +87,129 @@ public:
/// Get relative momentum error (no dimension).
double pt() const;
/// Get position in m
Vector_t R() const;
/// Get momentum
Vector_t P() const;
/// Get charge in Coulomb
double charge() const;
/// Get mass in GeV/c^2
double mass() const;
private:
// The particle's phase space coordinates:
double phase[6];
Vector_t R_m;
Vector_t P_m;
double charge_m;
double mass_m;
};
// Inline member functions.
// ------------------------------------------------------------------------
inline double &OpalParticle::operator[](int i) {
return phase[i];
inline double &OpalParticle::operator[](int i)
{
return i < PX? R_m[i] : P_m[i - PX];
}
inline double &OpalParticle::x()
{ return phase[X]; }
{
return R_m[X];
}
inline double &OpalParticle::y()
{ return phase[Y]; }
{
return R_m[Y];
}
inline double &OpalParticle::t()
{ return phase[T]; }
{
return R_m[T];
}
inline double &OpalParticle::px()
{ return phase[PX]; }
{
return P_m[X];
}
inline double &OpalParticle::py()
{ return phase[PY]; }
{
return P_m[Y];
}
inline double &OpalParticle::pt() {
return phase[PT];
inline double &OpalParticle::pt()
{
return P_m[T];
}
inline Vector_t &OpalParticle::R()
{
return R_m;
}
inline Vector_t &OpalParticle::P()
{
return P_m;
}
inline double OpalParticle::operator[](int i) const {
return phase[i];
inline double OpalParticle::operator[](int i) const
{
return i < PX? R_m[i]: P_m[i - PX];
}
inline double OpalParticle::x() const
{ return phase[X]; }
{
return R_m[X];
}
inline double OpalParticle::y() const
{ return phase[Y]; }
{
return R_m[Y];
}
inline double OpalParticle::t() const
{ return phase[T]; }
{
return R_m[T];
}
inline double OpalParticle::px() const
{ return phase[PX]; }
{
return P_m[X];
}
inline double OpalParticle::py() const
{ return phase[PY]; }
{
return P_m[Y];
}
inline double OpalParticle::pt() const
{
return P_m[T];
}
inline Vector_t OpalParticle::R() const
{
return R_m;
}
inline Vector_t OpalParticle::P() const
{
return P_m;
}
inline double OpalParticle::charge() const
{
return charge_m;
}
inline double OpalParticle::pt() const {
return phase[PT];
inline double OpalParticle::mass() const
{
return mass_m;
}
#endif // CLASSIC_OpalParticle_HH
#endif // CLASSIC_OpalParticle_HH
\ No newline at end of file
......@@ -93,7 +93,7 @@ void PartBunch::do_binaryRepart() {
get_bounds(rmin_m, rmax_m);
pbase_t* underlyingPbase =
dynamic_cast<pbase_t*>(pbase.get());
dynamic_cast<pbase_t*>(pbase_m.get());
BinaryRepartition(*underlyingPbase);
update();
......
......@@ -120,12 +120,12 @@ private:
//FIXME
ParticleLayout<double, 3> & getLayout() {
return pbase->getLayout();
return pbase_m->getLayout();
}
//FIXME
const ParticleLayout<double, 3>& getLayout() const {
return pbase->getLayout();
return pbase_m->getLayout();
}
};
......
......@@ -207,11 +207,11 @@ public:
Compatibility function push_back
*/
void push_back(OpalParticle p);
void push_back(OpalParticle const& p);
void set_part(FVector<double, 6> z, int ii);
void set_part(OpalParticle p, int ii);
void set_part(OpalParticle const& p, int ii);
OpalParticle get_part(int ii);
......@@ -441,11 +441,11 @@ public:
ParticleBConds<Position_t, Dimension>& getBConds() {
return pbase->getBConds();
return pbase_m->getBConds();
}
void setBConds(const ParticleBConds<Position_t, Dimension>& bc) {
pbase->setBConds(bc);
pbase_m->setBConds(bc);
}
bool singleInitNode() const;
......@@ -680,7 +680,7 @@ protected:
// flag to tell if we are a DC-beam
bool dcBeam_m;
double periodLength_m;
std::shared_ptr<AbstractParticle<T, Dim> > pbase;
std::shared_ptr<AbstractParticle<T, Dim> > pbase_m;
};
#include "PartBunchBase.hpp"
......
......@@ -92,7 +92,7 @@ PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartDat
dist_m(nullptr),
dcBeam_m(false),
periodLength_m(Physics::c / 1e9),
pbase(pb)
pbase_m(pb)
{
setup(pb);
}
......@@ -958,19 +958,14 @@ std::pair<Vector_t, double> PartBunchBase<T, Dim>::getLocalBoundingSphere() {
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::push_back(OpalParticle p) {
void PartBunchBase<T, Dim>::push_back(OpalParticle const& particle) {
Inform msg("PartBunch ");
size_t i = getLocalNum();
create