Commit ba000c54 authored by kraus's avatar kraus
Browse files

fixing striation bug

parent cf020342
......@@ -477,6 +477,41 @@ void Distribution::Create(size_t &numberOfParticles, double massIneV) {
INFOMSG("Distribution unknown." << endl;);
break;
}
if (emitting_m) {
unsigned int numAdditionalRNsPerParticle = 0;
if (emissionModel_m == EmissionModelT::ASTRA) {
numAdditionalRNsPerParticle = 2;
}
if (Options::cZero) {
numAdditionalRNsPerParticle *= 2;
}
int saveProcessor = -1;
for (size_t partIndex = 0; partIndex < numberOfParticles; ++ partIndex) {
// Save to each processor in turn.
++ saveProcessor;
if (saveProcessor >= Ippl::getNodes())
saveProcessor = 0;
if (Ippl::myNode() == saveProcessor) {
std::vector<double> rns;
for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
double x = gsl_rng_uniform(randGenEmit_m);
rns.push_back(x);
}
additionalRNs_m.push_back(rns);
} else {
for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
gsl_rng_uniform(randGenEmit_m);
}
}
}
}
// AAA Scale and shift coordinates according to distribution input.
ScaleDistCoordinates();
}
......@@ -927,7 +962,7 @@ void Distribution::AddDistributions() {
}
}
void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double &pz) {
void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
switch (emissionModel_m) {
......@@ -935,7 +970,7 @@ void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double
ApplyEmissModelNone(pz);
break;
case EmissionModelT::ASTRA:
ApplyEmissModelAstra(px, py, pz);
ApplyEmissModelAstra(px, py, pz, additionalRNs);
break;
case EmissionModelT::NONEQUIL:
ApplyEmissModelNonEquil(eZ, px, py, pz);
......@@ -945,10 +980,10 @@ void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double
}
}
void Distribution::ApplyEmissModelAstra(double &px, double &py, double &pz) {
void Distribution::ApplyEmissModelAstra(double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
double phi = 2.0 * acos(sqrt(gsl_rng_uniform(randGenEmit_m)));
double theta = 2.0 * Physics::pi * gsl_rng_uniform(randGenEmit_m);
double phi = 2.0 * acos(sqrt(additionalRNs[0]));
double theta = 2.0 * Physics::pi * additionalRNs[1];
px = pTotThermal_m * sin(phi) * cos(theta);
py = pTotThermal_m * sin(phi) * sin(theta);
......@@ -1253,7 +1288,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
- add flag in order to calculate tunes from FMLOWE to FMHIGHE
- eliminate physics and error
*/
std::string LineName = Attributes::getString(itsAttr[AttributesT::LINE]);
if (LineName != "") {
const BeamSequence* LineSequence = BeamSequence::find(LineName);
......@@ -1261,7 +1296,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
SpecificElementVisitor<Cyclotron> CyclotronVisitor(*LineSequence->fetchLine());
CyclotronVisitor.execute();
size_t NumberOfCyclotrons = CyclotronVisitor.size();
if (NumberOfCyclotrons > 1) {
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"I am confused, found more than one Cyclotron element in line");
......@@ -1271,7 +1306,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
"didn't find any Cyclotron element in line");
}
const Cyclotron* CyclotronElement = CyclotronVisitor.front();
*gmsg << "* ----------------------------------------------------" << endl;
*gmsg << "* About to find closed orbit and matched distribution " << endl;
*gmsg << "* I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
......@@ -1292,7 +1327,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
double lE,hE;
lE = fmLowE;
hE = fmHighE;
if ((lE<0) || (hE<0)) {
lE = E_m*1E-6;
hE = E_m*1E-6;
......@@ -1326,16 +1361,16 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[AttributesT::ORDERMAPS]),
writeMap);
if(siggen->match(Attributes::getReal(itsAttr[AttributesT::RESIDUUM]),
Attributes::getReal(itsAttr[AttributesT::MAXSTEPSSI]),
Attributes::getReal(itsAttr[AttributesT::MAXSTEPSCO]),
CyclotronElement->getPHIinit(),
Attributes::getReal(itsAttr[AttributesT::RGUESS]),
false)) {
std::array<double,3> Emit = siggen->getEmittances();
if (Attributes::getReal(itsAttr[AttributesT::RGUESS]) > 0)
*gmsg << "* RGUESS " << Attributes::getReal(itsAttr[AttributesT::RGUESS])/1000.0 << " (m) " << endl;
......@@ -1349,20 +1384,20 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
}
*gmsg << " \\\\" << endl;
}
/*
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: mm
moment: rad
*/
if(Options::cloTuneOnly)
throw OpalException("Do only CLO and tune calculation","... ");
auto sigma = siggen->getSigma();
// change units from mm to m
......@@ -1372,23 +1407,23 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
sigma(j, 2 * i) *= 1e-3;
}
}
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()", "Negative value on the diagonal of the sigma matrix.");
sigmaR_m[i] = std::sqrt(sigma(2 * i, 2 * i));
sigmaP_m[i] = std::sqrt(sigma(2 * i + 1, 2 * i + 1));
}
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
for (unsigned int i = 0; i < 3; ++ i) {
sigmaP_m[i] = ConverteVToBetaGamma(sigmaP_m[i], massIneV);
}
}
correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5)));
......@@ -1397,17 +1432,17 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
correlationMatrix_m(5, 0) = sigma(0, 5) / (sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (sqrt(sigma(1, 1) * sigma(5, 5)));
CreateDistributionGauss(numberOfParticles, massIneV);
CreateDistributionGauss(numberOfParticles, massIneV);
}
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
if (siggen)
delete siggen;
throw OpalException("Distribution::CreateMatchedGaussDistribution", "didn't find any matched distribution.");
}
if (siggen)
delete siggen;
}
......@@ -1658,6 +1693,16 @@ void Distribution::CreateOpalT(PartBunch &beam,
addedDistIt != addedDistributions_m.end(); addedDistIt++)
(*addedDistIt)->SetDistType();
// Check if this is to be an emitted beam.
CheckIfEmitted();
if (emitting_m) {
CheckEmissionParameters();
SetupEmissionModel(beam);
} else {
pmean_m = Vector_t(0.0, 0.0, ConverteVToBetaGamma(GetEkin(), beam.getM()));
}
/*
* If Options::cZero is true than we reflect generated distribution
* to ensure that the transverse averages are 0.0.
......@@ -1675,9 +1720,6 @@ void Distribution::CreateOpalT(PartBunch &beam,
*/
CalcPartPerDist(numberOfParticles);
// Check if this is to be an emitted beam.
CheckIfEmitted();
/*
* Force added distributions to the same emission condition as the main
* distribution.
......@@ -1708,13 +1750,6 @@ void Distribution::CreateOpalT(PartBunch &beam,
// Check number of particles in distribution.
CheckParticleNumber(numberOfParticles);
if (emitting_m) {
CheckEmissionParameters();
SetupEmissionModel(beam);
} else {
pmean_m = Vector_t(0.0, 0.0, ConverteVToBetaGamma(GetEkin(), beam.getM()));
}
/*
* Find max. and min. particle positions in the bunch. For the
* case of an emitted beam these will be in time (seconds). For
......@@ -1742,6 +1777,16 @@ void Distribution::CreateOpalT(PartBunch &beam,
InitializeBeam(beam);
WriteOutFileHeader();
if (emitting_m && Options::cZero) {
std::vector<std::vector<double> > mirrored;
const auto end = additionalRNs_m.end();
for (auto it = additionalRNs_m.begin(); it != end; ++ it) {
std::vector<double> tmp((*it).begin() + 2, (*it).end());
mirrored.push_back(tmp);
(*it).erase((*it).begin() + 2, (*it).end());
}
additionalRNs_m.insert(additionalRNs_m.end(), mirrored.begin(), mirrored.end());
}
/*
* If this is an injected beam, we create particles right away.
* Emitted beams get created during the course of the simulation.
......@@ -1832,7 +1877,11 @@ size_t Distribution::EmitParticles(PartBunch &beam, double eZ) {
double px = pxDist_m.at(particleIndex);
double py = pyDist_m.at(particleIndex);
double pz = pzDist_m.at(particleIndex);
ApplyEmissionModel(eZ, px, py, pz);
std::vector<double> additionalRNs;
if (additionalRNs_m.size() > particleIndex) {
additionalRNs = additionalRNs_m[particleIndex];
}
ApplyEmissionModel(eZ, px, py, pz, additionalRNs);
double particleGamma
= std::sqrt(1.0
......@@ -1886,6 +1935,11 @@ size_t Distribution::EmitParticles(PartBunch &beam, double eZ) {
std::swap(pyDist_m.at(*ptbErasedIt), pyDist_m.back());
std::swap(tOrZDist_m.at(*ptbErasedIt), tOrZDist_m.back());
std::swap(pzDist_m.at(*ptbErasedIt), pzDist_m.back());
if (additionalRNs_m.size() == xDist_m.size()) {
std::swap(additionalRNs_m.at(*ptbErasedIt), additionalRNs_m.back());
additionalRNs_m.pop_back();
}
xDist_m.pop_back();
pxDist_m.pop_back();
......@@ -4697,4 +4751,4 @@ void Distribution::WriteOutFileInjection() {
reduce(numberOfParticles, numberOfParticles, OpAddAssign());
}
}
}
}
\ No newline at end of file
......@@ -215,8 +215,8 @@ private:
// void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);
void AddDistributions();
void ApplyEmissionModel(double eZ, double &px, double &py, double &pz);
void ApplyEmissModelAstra(double &px, double &py, double &pz);
void ApplyEmissionModel(double eZ, double &px, double &py, double &pz, const std::vector<double> &additionalRNs);
void ApplyEmissModelAstra(double &px, double &py, double &pz, const std::vector<double> &additionalRNs);
void ApplyEmissModelNone(double &pz);
void ApplyEmissModelNonEquil(double eZ, double &px, double &py, double &pz);
void Create(size_t &numberOfParticles, double massIneV);
......@@ -332,6 +332,7 @@ private:
double cathodeTemp_m; /// Cathode temperature (K).
double emitEnergyUpperLimit_m; /// Upper limit on emission energy distribution (eV).
std::vector<std::vector<double> > additionalRNs_m;
// Beam coordinate containers.
std::vector<double> xDist_m;
......
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