Commit 8792c971 authored by kraus's avatar kraus
Browse files

adding more information to loss data sink

parent 93c63b24
......@@ -204,9 +204,11 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i],
dtime, intecoords, triId);
if(res >= 0) {
lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i],
itsBunch_m->ID[i], itsBunch_m->getT()*1e9,
turnnumber_m, itsBunch_m->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(itsBunch_m->ID[i],
itsBunch_m->R[i], itsBunch_m->P[i],
itsBunch_m->getT()*1e9,
itsBunch_m->Q[i], itsBunch_m->M[i]),
std::make_pair(turnnumber_m, itsBunch_m->bunchNum[i]));
itsBunch_m->Bin[i] = -1;
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
......@@ -589,7 +591,7 @@ void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
*gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;
elptr->initialise(itsBunch_m, elptr->getPScale());
BcParameter[1] = temperature;
BcParameter[2] = stop;
......
......@@ -75,7 +75,7 @@ bool CCollimator::doPreCheck(PartBunchBase<double, 3> *bunch) {
// rectangle collimators in cyclotron cylindrical coordinates
// when there is no particlematterinteraction, the particle hitting collimator is deleted directly
bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double /*t*/, const double /*tstep*/) {
bool flagNeedUpdate = false;
size_t tempnum = bunch->getLocalNum();
......@@ -88,7 +88,10 @@ bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumbe
/// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
if ((pflag != 0) && (bunch->Bin[i] != -1)) {
if (!parmatint_m)
lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
bunch->R[i], bunch->P[i],
bunch->getT(), bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
bunch->Bin[i] = -1;
flagNeedUpdate = true;
}
......
......@@ -408,8 +408,9 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &
}
if (flagNeedUpdate) {
lossDs_m->addParticle(RefPartBunch_m->R[id], RefPartBunch_m->P[id],
id, t, 0, RefPartBunch_m->bunchNum[id]);
lossDs_m->addParticle(OpalParticle(id, RefPartBunch_m->R[id], RefPartBunch_m->P[id],
t, RefPartBunch_m->Q[id], RefPartBunch_m->M[id]),
std::make_pair(0, RefPartBunch_m->bunchNum[id]));
RefPartBunch_m->Bin[id] = -1;
}
......@@ -438,23 +439,23 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/,
// Necessary for gap phase output -DW
if (0 <= tet && tet <= 45) waiting_for_gap = 1;
// dB_{z}/dr, dB_{z}/dtheta, B_{z}
double brint = 0.0, btint = 0.0, bzint = 0.0;
if ( this->interpolate(rad, tet_rad, brint, btint, bzint) ) {
/* Br */
double br = - brint * R[2];
/* Btheta */
double bt = - btint / rad * R[2];
/* Bz */
double bz = - bzint;
this->applyTrimCoil(rad, R[2], tet_rad, br, bz);
/* Br Btheta -> Bx By */
B[0] = br * std::cos(tet_rad) - bt * std::sin(tet_rad);
B[1] = br * std::sin(tet_rad) + bt * std::cos(tet_rad);
......@@ -697,11 +698,11 @@ bool Cyclotron::interpolate(const double& rad,
const double wr1 = xir - (double)ir;
// wr2 : the relative distance to the outer path radius
const double wr2 = 1.0 - wr1;
// the corresponding angle on the field map
// Note: this does not work if the start point of field map does not equal zero.
double tet_map = std::fmod(tet_rad * Physics::rad2deg, 360.0 / symmetry_m);
double xit = tet_map / BP.dtet;
int it = (int) xit;
......@@ -1578,4 +1579,4 @@ void Cyclotron::getFieldFromFile_Synchrocyclotron(const double &scaleFactor) {
void Cyclotron::getDimensions(double &/*zBegin*/, double &/*zEnd*/) const
{ }
#undef CHECK_CYC_FSCANF_EOF
#undef CHECK_CYC_FSCANF_EOF
\ No newline at end of file
......@@ -113,9 +113,11 @@ bool FlexibleCollimator::apply(const size_t &i, const double &t, Vector_t &/*E*/
if (pdead) {
if (lossDs_m) {
double frac = -R(2) / P(2) * recpgamma;
lossDs_m->addParticle(R, P,
RefPartBunch_m->ID[i],
t + frac * dt, 0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R, P,
t + frac * dt,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i], OpalParticle::SPATIAL),
std::make_pair(0, RefPartBunch_m->bunchNum[i]));
}
++ losses_m;
}
......@@ -291,4 +293,4 @@ void FlexibleCollimator::writeHolesAndQuadtree(const std::string &baseFilename)
out.close();
}
}
}
\ No newline at end of file
......@@ -81,8 +81,11 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t
dt * (R(2) + P(2) * recpgamma) > dt * middle) {
double frac = (middle - R(2)) / (P(2) * recpgamma);
lossDs_m->addParticle(R + frac * recpgamma * P,
P, RefPartBunch_m->ID[i], t + frac * dt, 0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R + frac * recpgamma * P, P,
t + frac * dt,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i],
OpalParticle::SPATIAL));
}
}
......@@ -117,11 +120,11 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
for (unsigned int i = 0; i < localNum; ++ i) {
const double recpgamma = Physics::c * dt / Util::getGamma(RefPartBunch_m->P[i]);
Vector_t shift = frac * recpgamma * RefPartBunch_m->P[i] - Vector_t(0, 0, middle);
lossDs_m->addParticle(RefPartBunch_m->R[i] + shift,
RefPartBunch_m->P[i],
RefPartBunch_m->ID[i],
time,
0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
RefPartBunch_m->R[i] + shift,
RefPartBunch_m->P[i],
time,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i]));
}
OpalData::OPENMODE openMode;
if (numPassages_m > 0) {
......
......@@ -108,8 +108,9 @@ bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const
// peak finder uses millimetre not metre
peakfinder_m->addParticle(probepoint * 1e3);
lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i], probepoint, bunch->P[i], t+dt,
bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
}
peakfinder_m->evaluate(turnnumber);
......
......@@ -106,7 +106,11 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
gmsgALL << getName() << ": particle " << id
<< " at " << refPartBunch_m->R[id]
<< " m out of the field map boundary" << endl;
lossDS_m->addParticle(refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id], id);
lossDS_m->addParticle(OpalParticle(id,
refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id],
t,
refPartBunch_m->Q[id], refPartBunch_m->M[id]),
std::make_pair(0, (short)0));
lossDS_m->save();
refPartBunch_m->Bin[id] = -1;
......
......@@ -78,7 +78,7 @@ bool Septum::doPreCheck(PartBunchBase<double, 3> *bunch) {
return false;
}
bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double /*t*/, const double /*tstep*/) {
bool flag = false;
const double slope = (yend_m - ystart_m) / (xend_m - xstart_m);
......@@ -100,7 +100,12 @@ bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/,
R(1) > ystart_m &&
R(1) < yend_m) {
lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
R, bunch->P[i],
bunch->getT(),
bunch->Q[i], bunch->M[i],
OpalParticle::SPATIAL),
std::make_pair(turnnumber, bunch->bunchNum[i]));
bunch->Bin[i] = -1;
flag = true;
}
......
......@@ -54,8 +54,12 @@ bool Source::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t &
if (online_m && R(2) <= 0.0 && P(2) < 0.0) {
double frac = -R(2) / (P(2) * recpgamma);
lossDs_m->addParticle(R + frac * recpgamma * P,
P, RefPartBunch_m->ID[i], t + frac * dt, 0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R + frac * recpgamma * P, P,
t + frac * dt,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i],
OpalParticle::SPATIAL),
std::make_pair(0, RefPartBunch_m->bunchNum[i]));
return true;
}
......
......@@ -131,8 +131,11 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
strippoint(0) = (B_m * B_m * bunch->R[i](0) - A_m * B_m* bunch->R[i](1) - A_m * C_m) / (R_m * R_m);
strippoint(1) = (A_m * A_m * bunch->R[i](1) - A_m * B_m* bunch->R[i](0) - B_m * C_m) / (R_m * R_m);
strippoint(2) = bunch->R[i](2);
lossDs_m->addParticle(strippoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
strippoint, bunch->P[i],
t+dt,
bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
flagNeedUpdate = true;
if (stop_m) {
......
......@@ -37,99 +37,31 @@ void DistributionMoments::compute(PartBunchBase<double, 3> const& bunch)
totalNumParticles_m = bunch.getTotalNum();
if (totalNumParticles_m == 0) {
return;
}
computeMoments(bunch);
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);
computeMoments(bunch.begin(), bunch.end());
}
void DistributionMoments::compute(std::vector<OpalParticle> const& particles)
void DistributionMoments::compute(const std::vector<OpalParticle>::const_iterator &first,
const std::vector<OpalParticle>::const_iterator &last)
{
reset();
unsigned int localNumParticles = particles.size();
unsigned int localNumParticles = last - first;
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);
computeMoments(first, last);
}
void DistributionMoments::computeMeanKineticEnergy(PartBunchBase<double, 3> const& bunch)
template<class InputIt>
void DistributionMoments::computeMoments(const InputIt &first, const InputIt &last)
{
double data[] = {0.0, 0.0};
for (OpalParticle const& particle: bunch) {
data[0] += Util::getKineticEnergy(particle.P(), particle.mass());
if (totalNumParticles_m == 0) {
return;
}
data[1] = bunch.getLocalNum();
allreduce(data, 2, std::plus<double>());
meanKineticEnergy_m = data[0] / data[1];
}
template<class Container>
void DistributionMoments::computeMoments(Container const& particles)
{
std::vector<double> localMoments(36);
for (OpalParticle const& particle: particles) {
for (InputIt it = first; it != last; ++ it) {
OpalParticle const& particle = *it;
unsigned int l = 6;
for (unsigned int i = 0; i < 6; ++ i) {
localMoments[i] += particle[i];
......@@ -143,8 +75,8 @@ void DistributionMoments::computeMoments(Container const& particles)
localMoments[l] += r2 * particle[i];
localMoments[l + 3] += r2 * r2;
}
double gamma = Util::getGamma(particle.P());
double eKin = (gamma - 1.0) * particle.mass();
double gamma = Util::getGamma(particle.getP());
double eKin = (gamma - 1.0) * particle.getMass();
localMoments[33] += eKin;
localMoments[34] += std::pow(eKin, 2);
localMoments[35] += gamma;
......@@ -152,11 +84,30 @@ void DistributionMoments::computeMoments(Container const& particles)
allreduce(localMoments.data(), localMoments.size(), std::plus<double>());
for (unsigned int i = 0; i < 6; ++ i) {
centroid_m[i] = localMoments[i];
fillMembers(localMoments);
}
void DistributionMoments::computeMeanKineticEnergy(PartBunchBase<double, 3> const& bunch)
{
double data[] = {0.0, 0.0};
for (OpalParticle const& particle: bunch) {
data[0] += Util::getKineticEnergy(particle.getP(), particle.getMass());
}
data[1] = bunch.getLocalNum();
allreduce(data, 2, std::plus<double>());
meanKineticEnergy_m = data[0] / data[1];
}
void DistributionMoments::fillMembers(std::vector<double> const& localMoments) {
Vector_t squaredEps, fac, squaredSumR, squaredSumP, sumRP;
double perParticle = 1.0 / totalNumParticles_m;
unsigned int l = 0;
for (; l < 6; ++ l) {
centroid_m[l] = localMoments[l];
}
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];
......@@ -164,13 +115,11 @@ void DistributionMoments::computeMoments(Container const& particles)
}
}
double perParticle = 1.0 / totalNumParticles_m;
for (unsigned int i = 0; i < 3; ++ i) {
for (unsigned int i = 0; i < 3; ++ i, ++ l) {
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 w3 = localMoments[l] * perParticle;
double w4 = localMoments[l + 3] * perParticle;
double tmp = w2 - std::pow(w1, 2);
halo_m(i) = (w4 + w1 * (-4 * w3 + 3 * w1 * (tmp + w2))) / tmp;
......@@ -180,6 +129,31 @@ void DistributionMoments::computeMoments(Container const& particles)
meanKineticEnergy_m = localMoments[33] * perParticle;
stdKineticEnergy_m = localMoments[34] * perParticle;
meanGamma_m = localMoments[35] * perParticle;
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::reset()
......
......@@ -15,7 +15,8 @@ class DistributionMoments {
public:
DistributionMoments();
void compute(std::vector<OpalParticle> const&);
void compute(const std::vector<OpalParticle>::const_iterator &,
const std::vector<OpalParticle>::const_iterator &);
void compute(PartBunchBase<double, 3> const&);
void computeMeanKineticEnergy(PartBunchBase<double, 3> const&);
......@@ -37,8 +38,9 @@ public:
private:
template<class Container>
void computeMoments(Container const&);
template<class InputIt>
void computeMoments(const InputIt &, const InputIt &);
void fillMembers(std::vector<double> const&);
void reset();
Vector_t meanR_m;
......
......@@ -22,12 +22,31 @@ OpalParticle::OpalParticle()
{}
OpalParticle::OpalParticle(double x, double px,
OpalParticle::OpalParticle(unsigned int id,
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),
double z, double pz,
double t,
double q, double m,
OpalParticle::Type type):
id_m(id),
R_m(x, y, z),
P_m(px, py, pz),
time_m(t),
charge_m(q),
mass_m(m)
mass_m(m),
type_m(type)
{}
OpalParticle::OpalParticle(unsigned int id,
Vector_t const& R, Vector_t const& P,
double t, double q, double m,
OpalParticle::Type type):
id_m(id),
R_m(R),
P_m(P),
time_m(t),
charge_m(q),
mass_m(m),
type_m(type)
{}
\ No newline at end of file
......@@ -25,191 +25,250 @@ class OpalParticle {
public:
// Particle coordinate numbers.
enum { X, Y, T, PX, PY, PT };
enum { X, Y, L, INVALID };
enum Type {
SPATIAL = 0,
TEMPORAL
};
/// Constructor.
// Construct particle with the given coordinates.
OpalParticle(double x, double px,
OpalParticle(unsigned int id,
double x, double px,
double y, double py,
double t, double pt,
double q, double m);
double z, double pz,
double time,
double q, double m, Type type = TEMPORAL);
OpalParticle(unsigned int id,
Vector_t const& R, Vector_t const& P,
double time, double q, double m, Type type = TEMPORAL);
OpalParticle();
/// Get coordinate.
// Access coordinate by index.
// Note above order of phase space coordinates!
double &operator[](int);
/// Set the horizontal position in m.
void setX(double) ;
/// Set the horizontal momentum.
void setPx(double);
/// Get reference to horizontal position in m.
double &x() ;
/// Set the vertical displacement in m.
void setY(double) ;
/// Get reference to horizontal momentum (no dimension).
double &px();
/// Set the vertical momentum.
void setPy(double);
/// Get reference to vertical displacement in m.
double &y() ;
/// Set longitudinal position in m.
void setZ(double) ;
/// Get reference to vertical momentum (no dimension).
double &py();
/// Set the longitudinal momentum
void setPz(double);
/// Get reference to longitudinal displacement c*t in m.
double &t() ;
/// Set position in m
void setR(Vector_t const&);
/// Get reference to relative momentum error (no dimension).
double &pt();
/// Set momentum