Commit 50e2a63e authored by kraus's avatar kraus
Browse files

add total charge and mass to the statistics

parent 8792c971
......@@ -26,6 +26,8 @@ DistributionMoments::DistributionMoments():
Dy_m(0.0),
DDy_m(0.0),
moments_m(),
totalCharge_m(0.0),
totalMass_m(0.0),
totalNumParticles_m(0)
{
std::fill(std::begin(centroid_m), std::end(centroid_m), 0.0);
......@@ -58,7 +60,7 @@ void DistributionMoments::computeMoments(const InputIt &first, const InputIt &la
return;
}
std::vector<double> localMoments(36);
std::vector<double> localMoments(38);
for (InputIt it = first; it != last; ++ it) {
OpalParticle const& particle = *it;
......@@ -80,6 +82,8 @@ void DistributionMoments::computeMoments(const InputIt &first, const InputIt &la
localMoments[33] += eKin;
localMoments[34] += std::pow(eKin, 2);
localMoments[35] += gamma;
localMoments[36] += particle.getCharge();
localMoments[37] += particle.getMass();
}
allreduce(localMoments.data(), localMoments.size(), std::plus<double>());
......@@ -129,6 +133,8 @@ void DistributionMoments::fillMembers(std::vector<double> const& localMoments) {
meanKineticEnergy_m = localMoments[33] * perParticle;
stdKineticEnergy_m = localMoments[34] * perParticle;
meanGamma_m = localMoments[35] * perParticle;
totalCharge_m = localMoments[36];
totalMass_m = localMoments[37];
for (unsigned int i = 0; i < 3; ++ i) {
meanR_m(i) = centroid_m[2 * i] * perParticle;
......
......@@ -35,7 +35,8 @@ public:
double getDDx() const;
double getDy() const;
double getDDy() const;
double getTotalCharge() const;
double getTotalMass() const;
private:
template<class InputIt>
......@@ -62,6 +63,8 @@ private:
double centroid_m[6];
FMatrix<double, 6, 6> moments_m;
double totalCharge_m;
double totalMass_m;
unsigned int totalNumParticles_m;
};
......@@ -155,4 +158,16 @@ double DistributionMoments::getDDy() const
return DDy_m;
}
inline
double DistributionMoments::getTotalCharge() const
{
return totalCharge_m;
}
inline
double DistributionMoments::getTotalMass() const
{
return totalMass_m;
}
#endif
\ No newline at end of file
......@@ -1228,138 +1228,10 @@ void PartBunchBase<T, Dim>::calcBeamParameters() {
IpplTimings::startTimer(statParamTimer_m);
momentsComputer_m.compute(*this);
// Vector_t eps2, fac, rsqsum, psqsum, rpsum;
// IpplTimings::startTimer(statParamTimer_m);
// const size_t locNp = getLocalNum();
// const size_t totalNum = getTotalNum();
// const double zero = 0.0;
// get_bounds(rmin_m, rmax_m);
// if(totalNum == 0) {
// for(unsigned int i = 0 ; i < Dim; i++) {
// rmean_m(i) = 0.0;
// pmean_m(i) = 0.0;
// rrms_m(i) = 0.0;
// prms_m(i) = 0.0;
// eps_norm_m(i) = 0.0;
// }
// rprms_m = 0.0;
// eKin_m = 0.0;
// eps_m = 0.0;
// IpplTimings::stopTimer(statParamTimer_m);
// return;
// }
// const size_t intN = calcMoments();
// const double N = static_cast<double>(intN);
// for(unsigned int i = 0 ; i < Dim; i++) {
// rmean_m(i) = centroid_m[2 * i] / N;
// pmean_m(i) = centroid_m[(2 * i) + 1] / N;
// rsqsum(i) = moments_m(2 * i, 2 * i) - N * rmean_m(i) * rmean_m(i);
// psqsum(i) = moments_m((2 * i) + 1, (2 * i) + 1) - N * pmean_m(i) * pmean_m(i);
// if(psqsum(i) < 0)
// psqsum(i) = 0;
// rpsum(i) = moments_m((2 * i), (2 * i) + 1) - N * rmean_m(i) * pmean_m(i);
// }
// eps2 = (rsqsum * psqsum - rpsum * rpsum) / (N * N);
// rpsum /= N;
// for(unsigned int i = 0 ; i < Dim; i++) {
// rrms_m(i) = std::sqrt(rsqsum(i) / N);
// prms_m(i) = std::sqrt(psqsum(i) / N);
// eps_norm_m(i) = std::sqrt(std::max(eps2(i), zero));
// double tmp = rrms_m(i) * prms_m(i);
// fac(i) = (tmp == 0) ? zero : 1.0 / tmp;
// }
// rprms_m = rpsum * fac;
// Dx_m = moments_m(0, 5) / N;
// DDx_m = moments_m(1, 5) / N;
// Dy_m = moments_m(2, 5) / N;
// DDy_m = moments_m(3, 5) / N;
// // Find unnormalized emittance.
// double gamma = 0.0;
// for(size_t i = 0; i < locNp; i++)
// gamma += std::sqrt(1.0 + dot(P[i], P[i]));
// allreduce(&gamma, 1, std::plus<double>());
// gamma /= N;
// calcEMean();
// // The computation of the energy spread is an estimation
// // based on the standard deviation of the longitudinal
// // momentum:
// // Var[f(P)] ~= (df/dP)(E[P])^2 Var[P]
// const double m0 = getM() * 1.E-6;
// double tmp = 1.0 / std::pow(eKin_m / m0 + 1., 2.0);
// if (OpalData::getInstance()->isInOPALCyclMode()) {
// dE_m = prms_m(1) * m0 * std::sqrt(1.0 - tmp);
// } else {
// dE_m = prms_m(2) * m0 * std::sqrt(1.0 - tmp);
// }
// eps_m = eps_norm_m / Vector_t(gamma * std::sqrt(1.0 - 1.0 / (gamma * gamma)));
IpplTimings::stopTimer(statParamTimer_m);
}
// template <class T, unsigned Dim>
// void PartBunchBase<T, Dim>::calcBeamParametersInitial() {
// const double N = static_cast<double>(getTotalNum());
// if(N == 0) {
// rmean_m = Vector_t(0.0);
// pmean_m = Vector_t(0.0);
// rrms_m = Vector_t(0.0);
// prms_m = Vector_t(0.0);
// eps_m = Vector_t(0.0);
// } else {
// if(Ippl::myNode() == 0) {
// // fixme: the following code is crap!
// // Only use one node as this function will get called only once before
// // particles have been emitted (at least in principle).
// Vector_t eps2, fac, rsqsum, psqsum, rpsum;
// const double zero = 0.0;
// const double N = static_cast<double>(pbin_m->getNp());
// calcMomentsInitial();
// for(unsigned int i = 0 ; i < Dim; i++) {
// rmean_m(i) = centroid_m[2 * i] / N;
// pmean_m(i) = centroid_m[(2 * i) + 1] / N;
// rsqsum(i) = moments_m(2 * i, 2 * i) - N * rmean_m(i) * rmean_m(i);
// psqsum(i) = moments_m((2 * i) + 1, (2 * i) + 1) - N * pmean_m(i) * pmean_m(i);
// if(psqsum(i) < 0)
// psqsum(i) = 0;
// rpsum(i) = moments_m((2 * i), (2 * i) + 1) - N * rmean_m(i) * pmean_m(i);
// }
// eps2 = (rsqsum * psqsum - rpsum * rpsum) / (N * N);
// rpsum /= N;
// for(unsigned int i = 0 ; i < Dim; i++) {
// rrms_m(i) = std::sqrt(rsqsum(i) / N);
// prms_m(i) = std::sqrt(psqsum(i) / N);
// eps_m(i) = std::sqrt(std::max(eps2(i), zero));
// double tmp = rrms_m(i) * prms_m(i);
// fac(i) = (tmp == 0) ? zero : 1.0 / tmp;
// }
// rprms_m = rpsum * fac;
// }
// }
// }
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getCouplingConstant() const {
return couplingConstant_m;
......
......@@ -203,24 +203,39 @@ void LossDataSink::writeHeaderH5() {
// Write file attributes to describe phase space to H5 file.
std::stringstream OPAL_version;
OPAL_version << OPAL_PROJECT_NAME << " " << OPAL_PROJECT_VERSION << " # git rev. " << Util::getGitRevision();
WRITE_FILEATTRIB_STRING ("OPAL_version", OPAL_version.str().c_str());
WRITE_FILEATTRIB_STRING ("tUnit", "s");
WRITE_FILEATTRIB_STRING ("xUnit", "m");
WRITE_FILEATTRIB_STRING ("yUnit", "m");
WRITE_FILEATTRIB_STRING ("zUnit", "m");
WRITE_FILEATTRIB_STRING ("pxUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING ("pyUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING ("pzUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING ("idUnit", "1");
WRITE_FILEATTRIB_STRING ("turnUnit", "1");
WRITE_FILEATTRIB_STRING ("timeUnit", "s");
WRITE_FILEATTRIB_STRING ("SPOSUnit", "mm");
WRITE_FILEATTRIB_STRING ("TIMEUnit", "s");
WRITE_FILEATTRIB_STRING ("mpart", "GeV");
WRITE_FILEATTRIB_STRING ("qi", "C");
WRITE_FILEATTRIB_STRING("OPAL_version", OPAL_version.str().c_str());
WRITE_FILEATTRIB_STRING("SPOSUnit", "m");
WRITE_FILEATTRIB_STRING("TIMEUnit", "s");
WRITE_FILEATTRIB_STRING("RefPartRUnit", "m");
WRITE_FILEATTRIB_STRING("RefPartPUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("GlobalTrackStepUnit", "1");
WRITE_FILEATTRIB_STRING("centroidUnit", "m");
WRITE_FILEATTRIB_STRING("RMSXUnit", "m");
WRITE_FILEATTRIB_STRING("MEANPUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("RMSPUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("#varepsilonUnit", "m rad");
WRITE_FILEATTRIB_STRING("#varepsilon-geomUnit", "m rad");
WRITE_FILEATTRIB_STRING("ENERGYUnit", "MeV");
WRITE_FILEATTRIB_STRING("dEUnit", "MeV");
WRITE_FILEATTRIB_STRING("TotalChargeUnit", "C");
WRITE_FILEATTRIB_STRING("TotalMassUnit", "MeV");
WRITE_FILEATTRIB_STRING("idUnit", "1");
WRITE_FILEATTRIB_STRING("xUnit", "m");
WRITE_FILEATTRIB_STRING("yUnit", "m");
WRITE_FILEATTRIB_STRING("zUnit", "m");
WRITE_FILEATTRIB_STRING("pxUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("pyUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("pzUnit", "#beta#gamma");
WRITE_FILEATTRIB_STRING("qUnit", "C");
WRITE_FILEATTRIB_STRING("mUnit", "MeV");
WRITE_FILEATTRIB_STRING("turnUnit", "1");
WRITE_FILEATTRIB_STRING("bunchNumUnit", "1");
WRITE_FILEATTRIB_STRING("timeUnit", "s");
}
void LossDataSink::addReferenceParticle(const Vector_t &x,
......@@ -356,23 +371,17 @@ void LossDataSink::saveH5(unsigned int setIdx) {
WRITE_STEPATTRIB_FLOAT64("#varepsilon-geom", (tmpVector = engine.getGeometricEmittance(), &tmpVector[0]), 3);
WRITE_STEPATTRIB_FLOAT64("ENERGY", (tmpDouble = engine.getMeanKineticEnergy(), &tmpDouble), 1);
WRITE_STEPATTRIB_FLOAT64("dE", (tmpDouble = engine.getStandardDeviationKineticEnergy(), &tmpDouble), 1);
double totalCharge = 0;
std::for_each(particles_m.begin() + startIdx,
particles_m.begin() + endIdx,
[&totalCharge](const OpalParticle& particle) { totalCharge += particle.getCharge(); });
WRITE_STEPATTRIB_FLOAT64("TotalCharge", &totalCharge, 1);
double totalMass = 0;
std::for_each(particles_m.begin() + startIdx,
particles_m.begin() + endIdx,
[&totalMass](const OpalParticle& particle) { totalMass += particle.getMass(); });
WRITE_STEPATTRIB_FLOAT64("TotalMass", &totalMass, 1);
WRITE_STEPATTRIB_FLOAT64("TotalCharge", (tmpDouble = engine.getTotalCharge(), &tmpDouble), 1);
WRITE_STEPATTRIB_FLOAT64("TotalMass", (tmpDouble = engine.getTotalMass(), &tmpDouble), 1);
// Write all data
std::vector<char> buffer(nLoc * sizeof(h5_float64_t));
h5_float64_t *f64buffer = reinterpret_cast<h5_float64_t*>(&buffer[0]);
h5_int64_t *i64buffer = reinterpret_cast<h5_int64_t*>(&buffer[0]);
::i64transform(particles_m, startIdx, nLoc, i64buffer,
[](const OpalParticle &particle) { return particle.getId(); });
WRITE_DATA_INT64 ("id", i64buffer);
::f64transform(particles_m, startIdx, nLoc, f64buffer,
[](const OpalParticle &particle) { return particle.getX(); });
WRITE_DATA_FLOAT64 ("x", f64buffer);
......@@ -391,9 +400,6 @@ void LossDataSink::saveH5(unsigned int setIdx) {
::f64transform(particles_m, startIdx, nLoc, f64buffer,
[](const OpalParticle &particle) { return particle.getPz(); });
WRITE_DATA_FLOAT64 ("pz", f64buffer);
::i64transform(particles_m, startIdx, nLoc, i64buffer,
[](const OpalParticle &particle) { return particle.getId(); });
WRITE_DATA_INT64 ("id", i64buffer);
::f64transform(particles_m, startIdx, nLoc, f64buffer,
[](const OpalParticle &particle) { return particle.getCharge(); });
WRITE_DATA_FLOAT64 ("q", f64buffer);
......
......@@ -194,7 +194,6 @@ void H5PartWrapperForPT::writeHeader() {
OPAL_version << OPAL_PROJECT_NAME << " " << OPAL_PROJECT_VERSION << " # git rev. " << Util::getGitRevision();
WRITESTRINGFILEATTRIB(file_m, "OPAL_version", OPAL_version.str().c_str());
WRITESTRINGFILEATTRIB(file_m, "idUnit", "1");
WRITESTRINGFILEATTRIB(file_m, "xUnit", "m");
WRITESTRINGFILEATTRIB(file_m, "yUnit", "m");
......@@ -499,4 +498,4 @@ void H5PartWrapperForPT::writeStepData(PartBunchBase<double, 3>* bunch) {
reportOnError(herr, __FILE__, __LINE__);
}
}
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment