Commit 79d7c0c6 authored by adelmann's avatar adelmann 🎗
Browse files

correct Ekin to match energy in input file

parent d5058870
...@@ -2067,6 +2067,23 @@ void PartBunch::calcBeamParameters_cycl() { ...@@ -2067,6 +2067,23 @@ void PartBunch::calcBeamParameters_cycl() {
eps_m = eps_norm_m / Vector_t(betagamma); eps_m = eps_norm_m / Vector_t(betagamma);
} }
void PartBunch::correctEnergy(double avrgp_m) {
const double totalNp = static_cast<double>(this->getTotalNum());
const double locNp = static_cast<double>(this->getLocalNum());
double avrgp = 0.0;
for(unsigned int k = 0; k < locNp; k++)
avrgp += sqrt(dot(P[k], P[k]));
reduce(avrgp, avrgp, OpAddAssign());
avrgp /= totalNp;
for(unsigned int k = 0; k < locNp; k++)
P[k](2) = P[k](2) - avrgp + avrgp_m;
}
void PartBunch::calcEMean() { void PartBunch::calcEMean() {
const double totalNp = static_cast<double>(this->getTotalNum()); const double totalNp = static_cast<double>(this->getTotalNum());
...@@ -2080,17 +2097,6 @@ void PartBunch::calcEMean() { ...@@ -2080,17 +2097,6 @@ void PartBunch::calcEMean() {
reduce(eKin_m, eKin_m, OpAddAssign()); reduce(eKin_m, eKin_m, OpAddAssign());
eKin_m /= totalNp; eKin_m /= totalNp;
/*
std::for_each(P.begin(), P.end(), [&](Vector_t x) { pm += x; });
double bega = std::sqrt(dot(pm,pm));
reduce(bega, bega, OpAddAssign());
bega /= totalNp;
double ga = std::sqrt(bega*bega + 1.0);
eKin_m = (ga-1.0) * getM() * 1e-6;
*/
} }
......
...@@ -517,10 +517,10 @@ private: ...@@ -517,10 +517,10 @@ private:
double calculateAngle(double x, double y); double calculateAngle(double x, double y);
double calculateAngle2(double x, double y); double calculateAngle2(double x, double y);
public:
double pSqr(Vector_t p1, Vector_t p2) const;
void calcEMean(); // update eKin_m; void calcEMean(); // update eKin_m;
void correctEnergy(double avrgp);
private:
/* /*
Member variables starts here Member variables starts here
...@@ -682,12 +682,6 @@ private: ...@@ -682,12 +682,6 @@ private:
}; };
inline double PartBunch::pSqr(Vector_t p1, Vector_t p2) const {
return p1[0]*p2[0] + p1[0]*p2[1] + p1[0]*p2[2] + p1[1]*p2[0] + p1[1]*p2[1] + p1[1]*p2[2] + p1[2]*p2[0] + p1[2]*p2[1] + p1[2]*p2[2];
}
inline inline
bool PartBunch::isDcBeam() { return dcBeam_m;} bool PartBunch::isDcBeam() { return dcBeam_m;}
......
...@@ -1628,7 +1628,7 @@ void Distribution::CreateDistributionGauss(size_t numberOfParticles, double mass ...@@ -1628,7 +1628,7 @@ void Distribution::CreateDistributionGauss(size_t numberOfParticles, double mass
GenerateTransverseGauss(numberOfParticles); GenerateTransverseGauss(numberOfParticles);
GenerateLongFlattopT(numberOfParticles); GenerateLongFlattopT(numberOfParticles);
} else { } else {
GenerateGaussZ(numberOfParticles); GenerateGaussZ(numberOfParticles);
} }
} }
...@@ -2714,7 +2714,6 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) { ...@@ -2714,7 +2714,6 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
gsl_matrix_set (m, i, j, 0.0); gsl_matrix_set (m, i, j, 0.0);
} }
} }
gsl_matrix_set (m,0,1, distCorr_m.at(0)); gsl_matrix_set (m,0,1, distCorr_m.at(0));
gsl_matrix_set (m,1,0, distCorr_m.at(0)); gsl_matrix_set (m,1,0, distCorr_m.at(0));
gsl_matrix_set (m,2,3, distCorr_m.at(1)); gsl_matrix_set (m,2,3, distCorr_m.at(1));
...@@ -2726,100 +2725,100 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) { ...@@ -2726,100 +2725,100 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
gsl_matrix_set (m,0,5, distCorr_m.at(3)); gsl_matrix_set (m,0,5, distCorr_m.at(3));
gsl_matrix_set (m,5,1, distCorr_m.at(4)); gsl_matrix_set (m,5,1, distCorr_m.at(4));
gsl_matrix_set (m,1,5, distCorr_m.at(4)); gsl_matrix_set (m,1,5, distCorr_m.at(4));
gsl_matrix_set (m,4,0, distCorr_m.at(5)); gsl_matrix_set (m,4,0, distCorr_m.at(5));
gsl_matrix_set (m,0,4, distCorr_m.at(5)); gsl_matrix_set (m,0,4, distCorr_m.at(5));
gsl_matrix_set (m,4,1, distCorr_m.at(6)); gsl_matrix_set (m,4,1, distCorr_m.at(6));
gsl_matrix_set (m,1,4, distCorr_m.at(6)); gsl_matrix_set (m,1,4, distCorr_m.at(6));
#define DISTDBG1
#ifdef DISTDBG1 #ifdef DISTDBG1
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
for (int j = 0; j < 6; j++) { for (int j = 0; j < 6; j++) {
if (j==0) if (j==0)
*gmsg << "* r(" << std::setprecision(1) << i << "," << std::setprecision(1) << j << ") = " *gmsg << "* r(" << std::setprecision(1) << i << "," << std::setprecision(1) << j << ") = "
<< std::setprecision(3) << gsl_matrix_get (m, i, j); << std::setprecision(3) << gsl_matrix_get (m, i, j);
else else
*gmsg << "\t" << std::setprecision(3) << gsl_matrix_get (m, i, j); *gmsg << "\t" << std::setprecision(3) << gsl_matrix_get (m, i, j);
} }
*gmsg << endl; *gmsg << endl;
} }
#endif #endif
int errcode = gsl_linalg_cholesky_decomp(m); int errcode = gsl_linalg_cholesky_decomp(m);
if (errcode == GSL_EDOM) { if (errcode == GSL_EDOM) {
INFOMSG("gsl_linalg_cholesky_decomp faliled" << endl); INFOMSG("gsl_linalg_cholesky_decomp faliled" << endl);
exit(1); exit(1);
} }
// so make the lower part zero. // so make the lower part zero.
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
for (int j = 0; j < i ; j++) { for (int j = 0; j < i ; j++) {
gsl_matrix_set (m, i, j, 0.0); gsl_matrix_set (m, i, j, 0.0);
} }
} }
gsl_matrix_transpose(m); gsl_matrix_transpose(m);
#define DISTDBG2
#ifdef DISTDBG2 #ifdef DISTDBG2
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
for (int j = 0; j < 6; j++) { for (int j = 0; j < 6; j++) {
if (j==0) if (j==0)
*gmsg << "* r(" << std::setprecision(1) << i << "," << std::setprecision(1) << j << ") = " *gmsg << "* r(" << std::setprecision(1) << i << "," << std::setprecision(1) << j << ") = "
<< std::setprecision(3) << gsl_matrix_get (m, i, j); << std::setprecision(3) << gsl_matrix_get (m, i, j);
else else
*gmsg << "\t" << std::setprecision(3) << gsl_matrix_get (m, i, j); *gmsg << "\t" << std::setprecision(3) << gsl_matrix_get (m, i, j);
} }
*gmsg << endl; *gmsg << endl;
} }
#endif #endif
int saveProcessor = -1; int saveProcessor = -1;
for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) { for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
double rval = 0.0; double rval = 0.0;
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[0]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[0])
; ;
gsl_vector_set(rx,0,rval); gsl_vector_set(rx,0,rval);
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[0]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[0])
; ;
gsl_vector_set(rx,1,rval); gsl_vector_set(rx,1,rval);
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[1]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[1])
; ;
gsl_vector_set(rx,2,rval); gsl_vector_set(rx,2,rval);
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[1]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[1])
; ;
gsl_vector_set(rx,3,rval); gsl_vector_set(rx,3,rval);
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[2]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffR_m[2])
; ;
gsl_vector_set(rx,4,rval); gsl_vector_set(rx,4,rval);
while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[2]) while (std::abs((rval = gsl_ran_gaussian (randGen,1.0)))>cutoffP_m[2])
; ;
gsl_vector_set(rx,5,rval); gsl_vector_set(rx,5,rval);
// Save to each processor in turn. // Save to each processor in turn.
saveProcessor++; saveProcessor++;
if (saveProcessor >= Ippl::getNodes()) if (saveProcessor >= Ippl::getNodes())
saveProcessor = 0; saveProcessor = 0;
if (Ippl::myNode() == saveProcessor) { if (Ippl::myNode() == saveProcessor) {
if (gsl_blas_dgemv(CblasNoTrans,1.0,m,rx,0.0,ry)) { if (gsl_blas_dgemv(CblasNoTrans,1.0,m,rx,0.0,ry)) {
INFOMSG("oops... something wrong with GSL matvec\n"); INFOMSG("oops... something wrong with GSL matvec\n");
exit(1); exit(1);
} }
xDist_m.push_back( sigmaR_m[0]*gsl_vector_get(ry, 0)); xDist_m.push_back( sigmaR_m[0]*gsl_vector_get(ry, 0));
pxDist_m.push_back( sigmaP_m[0]*gsl_vector_get(ry, 1)); pxDist_m.push_back( sigmaP_m[0]*gsl_vector_get(ry, 1));
yDist_m.push_back( sigmaR_m[1]*gsl_vector_get(ry, 2)); yDist_m.push_back( sigmaR_m[1]*gsl_vector_get(ry, 2));
pyDist_m.push_back( sigmaP_m[1]*gsl_vector_get(ry, 3)); pyDist_m.push_back( sigmaP_m[1]*gsl_vector_get(ry, 3));
tOrZDist_m.push_back( sigmaR_m[2]*gsl_vector_get(ry, 4)); tOrZDist_m.push_back( sigmaR_m[2]*gsl_vector_get(ry, 4));
pzDist_m.push_back(avrgpz_m +(sigmaP_m[2]*gsl_vector_get(ry, 5))); pzDist_m.push_back(avrgpz_m +(sigmaP_m[2]*gsl_vector_get(ry, 5)));
} }
} }
/* /*
//std::for_each(v.rbegin(), v.rend(), [&](int n) { sum_of_elements += n; }); //std::for_each(v.rbegin(), v.rend(), [&](int n) { sum_of_elements += n; });
double pxm = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0); double pxm = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0);
...@@ -3133,7 +3132,9 @@ void Distribution::InjectBeam(PartBunch &beam) { ...@@ -3133,7 +3132,9 @@ void Distribution::InjectBeam(PartBunch &beam) {
beam.boundp(); beam.boundp();
beam.correctEnergy(avrgpz_m);
beam.calcEMean();
} }
bool Distribution::GetIfDistEmitting() { bool Distribution::GetIfDistEmitting() {
......
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