Commit 53390248 authored by Daniel Winklehner's avatar Daniel Winklehner
Browse files

1. Changed the way eKin_m is calculated in PartBunch to be consistent between...

1. Changed the way eKin_m is calculated in PartBunch to be consistent between P and E especially when restarting, 2. Removed some debug outputs.
parent 1801d146
......@@ -2307,17 +2307,27 @@ void PartBunch::correctEnergy(double avrgp_m) {
void PartBunch::calcEMean() {
const double totalNp = static_cast<double>(this->getTotalNum());
const double locNp = static_cast<double>(this->getLocalNum());
const double totalNp = static_cast<double>(this->getTotalNum());
const double locNp = static_cast<double>(this->getLocalNum());
Vector_t pm(0.0);
Vector_t meanP_temp = Vector_t(0.0);
eKin_m = 0.0;
for(unsigned int k = 0; k < locNp; k++)
eKin_m += (sqrt(dot(P[k], P[k]) + 1.0) - 1.0) * getM() * 1e-6;
reduce(eKin_m, eKin_m, OpAddAssign());
eKin_m /= totalNp;
eKin_m = 0.0;
for(unsigned int k = 0; k < locNp; k++)
meanP_temp += P[k];
//eKin_m += sqrt(dot(P[k], P[k]) + 1.0);
//eKin_m -= locNp;
//eKin_m *= getM() * 1.0e-6;
reduce(meanP_temp, meanP_temp, OpAddAssign());
//reduce(eKin_m, eKin_m, OpAddAssign());
meanP_temp /= totalNp;
//eKin_m /= totalNp;
eKin_m = (sqrt(dot(meanP_temp, meanP_temp) + 1.0) - 1.0) * getM() * 1.0e-6;
}
......
......@@ -376,6 +376,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
// Note: Nothing else has to be set, b/c everything comes from the h5 file -DW
}
/*
// TEMP Debug Output -DW
Vector_t const meanP = calcMeanP();
*gmsg << endl;
......@@ -387,6 +388,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << "referencePt = " << referencePt << endl;
*gmsg << endl;
// ENDTEMP
*/
sinRefTheta_m = sin(referenceTheta / 180.0 * pi);
cosRefTheta_m = cos(referenceTheta / 180.0 * pi);
......@@ -400,7 +402,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << "* Bunch global starting momenta:" << endl;
*gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
*gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
*gmsg << "* Total reference momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
......@@ -3457,6 +3459,7 @@ Vector_t ParallelCyclotronTracker::calcMeanR() const {
Vector_t ParallelCyclotronTracker::calcMeanP() const {
Vector_t meanP(0.0, 0.0, 0.0);
for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
for(int d = 0; d < 3; ++d) {
meanP(d) += itsBunch->P[i](d);
......@@ -4150,6 +4153,50 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceStatData() {
double const betagamma_temp = sqrt(dot(meanP, meanP));
// TEMP DEBUG output -DW
//*gmsg << endl;
//*gmsg << "* BetaGamma from meanP = " << betagamma_temp << endl;
//*gmsg << "* BetaGamma from meanE = " << sqrt((1.0 + E / itsBunch->getM() * 1.0e6) * (1.0 + E / itsBunch->getM() * 1.0e6) - 1.0) << endl;
/*
*gmsg << endl;
*gmsg << "* Comparison of beta gamma from meanE and meanP respectively:" << endl;
const double totalNp = itsBunch->getTotalNum();
const double locNp = itsBunch->getLocalNum();
Vector_t meanP_temp = Vector_t(0.0, 0.0, 0.0);
double meanE_temp = 0.0;
for(unsigned int k = 0; k < locNp; ++k) {
meanE_temp += sqrt(dot(itsBunch->P[k], itsBunch->P[k]) + 1.0);
for(int d = 0; d < 3; ++d) {
meanP_temp(d) += itsBunch->P[k](d);
}
}
meanE_temp -= locNp;
meanE_temp *= itsBunch->getM() * 1.0e-6;
reduce(meanE_temp, meanE_temp, OpAddAssign());
reduce(meanP_temp, meanP_temp, OpAddAssign());
meanP_temp /= totalNp;
meanE_temp /= totalNp;
*gmsg << "* meanE from meanE = " << meanE_temp << endl;
*gmsg << "* meanE from meanP = " << (sqrt(dot(meanP_temp, meanP_temp) + 1.0) - 1.0) * itsBunch->getM() * 1.0e-6 << endl;
*gmsg << "* BetaGamma from meanP = " << sqrt(dot(meanP_temp, meanP_temp)) << endl;
*gmsg << "* BetaGamma from meanE = " << sqrt((1.0 + meanE_temp / itsBunch->getM() * 1.0e6) * (1.0 + meanE_temp / itsBunch->getM() * 1.0e6) - 1.0) << endl;
// END TEMP
*/
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
double const theta = atan2(meanR(1), meanR(0));
......
......@@ -786,15 +786,13 @@ void Distribution::DoRestartOpalCycl(PartBunch &beam, size_t Np, int restartStep
*gmsg << "* Restart Energy = " << meanE << " (MeV), Path lenght = " << pathLength << " (m)" << endl;
*gmsg << "beam.getM() = " << beam.getM() << endl;
// beam.getM() in eV, meanE in MeV had to change 1e3 to 1e6 -DW
double ga = 1 + meanE/beam.getM()*1.0e6;
double be = sqrt(1.0-(1.0/(ga*ga)));
double ga = 1 + meanE / beam.getM() * 1.0e6;
double be = sqrt(1.0 - (1.0 / (ga * ga)));
bega_m = be*ga;
bega_m = be * ga;
*gmsg << "* Restart Energy = " << meanE << " (MeV), gamma = " << ga << ", beta = " << be << endl;
*gmsg << "* Gamma = " << ga << ", Beta = " << be << endl;
std::unique_ptr<char[]> varray(new char[(localN)*sizeof(double)]);
......
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