diff --git a/src/Classic/AbsBeamline/ElementBase.h b/src/Classic/AbsBeamline/ElementBase.h
index 86affdbaaafc72ad6ea57eaab0782e27a4a79d5e..181e58d701a2db904d2b8715bed8cd767e57e369 100644
--- a/src/Classic/AbsBeamline/ElementBase.h
+++ b/src/Classic/AbsBeamline/ElementBase.h
@@ -222,6 +222,12 @@ public:
// This may be the arc length or the straight length.
virtual void setElementLength(double length);
+ virtual void getElementDimensions(double &begin,
+ double &end) const {
+ begin = 0.0;
+ end = getElementLength();
+ }
+
/// Get origin position.
// Return the arc length from the entrance to the origin of the element
// (origin >= 0)
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index 59ddca0e01892491e75b8cbbab5679eb3a16d1e3..1ec1b01426578f6b241066572c3ac689398a114a 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -171,6 +171,9 @@ void RFCavity::accept(BeamlineVisitor &visitor) const {
void RFCavity::addKR(int i, double t, Vector_t &K) {
double pz = RefPartBunch_m->getZ(i) - startField_m;
+ if (pz < 0.0 ||
+ pz >= length_m) return;
+
const Vector_t tmpR(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);
double k = -Physics::q_e / (2.0 * RefPartBunch_m->getGamma(i) * Physics::EMASS);
@@ -183,6 +186,9 @@ void RFCavity::addKR(int i, double t, Vector_t &K) {
double wtf = frequency_m * t + phase_m;
double kj = k * scale_m * (tmpE(2) * cos(wtf) - RefPartBunch_m->getBeta(i) * frequency_m * Ez * sin(wtf) / Physics::c);
+
+ *gmsg << __DBGMSG__ << std::scientific << tmpE(2) << "\t" << Ez << "\t" << kj << "\t" << frequency_m << "\t" << scale_m << "\t" << k << endl;
+
K += Vector_t(kj, kj, 0.0);
}
@@ -196,12 +202,15 @@ void RFCavity::addKT(int i, double t, Vector_t &K) {
RefPartBunch_m->actT();
+ double pz = RefPartBunch_m->getZ(i) - startField_m;
+ if (pz < 0.0 ||
+ pz >= length_m) return;
+
//XXX: BET parameter, default is 1.
//If cxy != 1, then cxy = true
bool cxy = false; // default
double kx = 0.0, ky = 0.0;
if(cxy) {
- double pz = RefPartBunch_m->getZ(i) - startField_m;
const Vector_t tmpA(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
@@ -834,4 +843,21 @@ bool RFCavity::isInside(const Vector_t &r) const {
}
return false;
+}
+
+double RFCavity::getElementLength() const {
+ double length = ElementBase::getElementLength();
+ if (length < 1e-10 && fieldmap_m != NULL) {
+ double start, end, tmp;
+ fieldmap_m->getFieldDimensions(start, end, tmp, tmp);
+ length = end - start;
+ }
+
+ return length;
+}
+
+void RFCavity::getElementDimensions(double &begin,
+ double &end) const {
+ double tmp;
+ fieldmap_m->getFieldDimensions(begin, end, tmp, tmp);
}
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/RFCavity.h b/src/Classic/AbsBeamline/RFCavity.h
index 14cd5f0c264002cfa5c5302b899ac8c933035547..207817a828dacb6257a692affbb08af862a2c008 100644
--- a/src/Classic/AbsBeamline/RFCavity.h
+++ b/src/Classic/AbsBeamline/RFCavity.h
@@ -203,6 +203,11 @@ public:
void setFrequencyModelName(std::string name);
std::string getFrequencyModelName();
+ virtual double getElementLength() const;
+ virtual void getElementDimensions(double &begin,
+ double &end) const;
+
+
protected:
std::shared_ptr<AbstractTimeDependence> phase_td_m;
std::string phase_name_m;
@@ -505,4 +510,4 @@ std::string RFCavity::getFrequencyModelName() {
return frequency_name_m;
}
-#endif // CLASSIC_RFCavity_HH
+#endif // CLASSIC_RFCavity_HH
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Solenoid.cpp b/src/Classic/AbsBeamline/Solenoid.cpp
index 011ccf5b2e1808f0fcfaa01c5ee8f9aa9e37e1f8..ed8387f94b9e11e3e31269b9888c8804553d513e 100644
--- a/src/Classic/AbsBeamline/Solenoid.cpp
+++ b/src/Classic/AbsBeamline/Solenoid.cpp
@@ -101,9 +101,12 @@ bool Solenoid::getFast() const {
void Solenoid::addKR(int i, double t, Vector_t &K) {
Inform msg("Solenoid::addKR()");
+ double pz = RefPartBunch_m->getZ(i) - startField_m;
+ if (pz < 0.0 ||
+ pz >= length_m) return;
+
Vector_t tmpE(0.0, 0.0, 0.0);
Vector_t tmpB(0.0, 0.0, 0.0);
- double pz = RefPartBunch_m->getZ(i) - startField_m;
const Vector_t tmpA(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);
myFieldmap_m->getFieldstrength(tmpA, tmpE, tmpB);
@@ -121,11 +124,14 @@ void Solenoid::addKT(int i, double t, Vector_t &K) {
Inform msg("Solenoid::addKT()");
double dbdz, emg;
+ double pz = RefPartBunch_m->getZ(i) - startField_m;
+ if (pz < 0.0 ||
+ pz >= length_m) return;
+
Vector_t tmpE(0.0, 0.0, 0.0);
Vector_t tmpB(0.0, 0.0, 0.0);
Vector_t tmpE_diff(0.0, 0.0, 0.0);
Vector_t tmpB_diff(0.0, 0.0, 0.0);
- double pz = RefPartBunch_m->getZ(i) - startField_m;
const Vector_t tmpA(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);
// define z direction:
@@ -253,3 +259,14 @@ ElementBase::ElementType Solenoid::getType() const {
bool Solenoid::isInside(const Vector_t &r) const {
return (r(2) >= startField_m && r(2) < startField_m + length_m && isInsideTransverse(r) && myFieldmap_m->isInside(r));
}
+
+
+double Solenoid::getElementLength() const {
+ return length_m;
+}
+
+void Solenoid::getElementDimensions(double &begin,
+ double &end) const {
+ begin = startField_m;
+ end = begin + length_m;
+}
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Solenoid.h b/src/Classic/AbsBeamline/Solenoid.h
index 705b4e6175a41d34e6b3c91538aff70cbff6bbb7..74565c97e517f42b3cdba03ebf036242ef2f0f0c 100644
--- a/src/Classic/AbsBeamline/Solenoid.h
+++ b/src/Classic/AbsBeamline/Solenoid.h
@@ -85,6 +85,10 @@ public:
virtual void getDimensions(double &zBegin, double &zEnd) const;
virtual bool isInside(const Vector_t &r) const;
+
+ virtual double getElementLength() const;
+
+ virtual void getElementDimensions(double &zBegin, double &zEnd) const;
private:
// std::string name; /**< The name of the object*/
@@ -101,4 +105,4 @@ private:
void operator=(const Solenoid &);
};
-#endif // CLASSIC_Solenoid_HH
+#endif // CLASSIC_Solenoid_HH
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/TravelingWave.cpp b/src/Classic/AbsBeamline/TravelingWave.cpp
index 98dbef1e7c7fd78c9d285dd85260ce0c2f36637e..5e381656877b5e79be1dbf18174b3e7549787488 100644
--- a/src/Classic/AbsBeamline/TravelingWave.cpp
+++ b/src/Classic/AbsBeamline/TravelingWave.cpp
@@ -427,6 +427,16 @@ void TravelingWave::getDimensions(double &zBegin, double &zEnd) const {
}
+double TravelingWave::getElementLength() const {
+ return length_m;
+}
+
+void TravelingWave::getElementDimensions(double &begin,
+ double &end) const {
+ begin = -0.5 * PeriodLength_m;
+ end = begin + length_m;
+}
+
ElementBase::ElementType TravelingWave::getType() const {
return TRAVELINGWAVE;
}
diff --git a/src/Classic/AbsBeamline/TravelingWave.h b/src/Classic/AbsBeamline/TravelingWave.h
index 7e6dd94710b6d881e6ad5a9e9e88cf671ad1d010..a349aa3bd81237d4f1d6530c3475bc99cfa7b8d7 100644
--- a/src/Classic/AbsBeamline/TravelingWave.h
+++ b/src/Classic/AbsBeamline/TravelingWave.h
@@ -98,6 +98,12 @@ public:
virtual void getDimensions(double &zBegin, double &zEnd) const;
virtual bool isInside(const Vector_t &r) const;
+
+ virtual double getElementLength() const;
+ virtual void getElementDimensions(double &begin,
+ double &end) const;
+
+
private:
Fieldmap *CoreFieldmap_m;
/* Fieldmap *EntryFringeField_m; */
@@ -231,4 +237,4 @@ inline
void TravelingWave::setMode(double mode) {
Mode_m = mode;
}
-#endif // CLASSIC_TravelingWave_HH
+#endif // CLASSIC_TravelingWave_HH
\ No newline at end of file
diff --git a/src/Classic/Fields/FM1DProfile1.cpp b/src/Classic/Fields/FM1DProfile1.cpp
index d607cf069b0978882a06f714484ab417987c0071..2862cb330b479b83aeeb0f08f4f584c463c0c591 100644
--- a/src/Classic/Fields/FM1DProfile1.cpp
+++ b/src/Classic/Fields/FM1DProfile1.cpp
@@ -190,17 +190,19 @@ void FM1DProfile1::readMap() {
}
- if (computeEntranceFringe(entranceParameter1_m) < computeEntranceFringe(entranceParameter3_m))
- throw GeneralClassicException("FM1DProfile1::readMap",
- "The entry fringe field should be defined such that\n"
- "the field is bigger at z = 'Entrance Parameter 1' than at\n"
- "z = 'Entrance Parameter 3'");
-
- if (computeExitFringe(exitParameter1_m) < computeExitFringe(exitParameter3_m))
- throw GeneralClassicException("FM1DProfile1::readMap",
- "The exit fringe field should be defined such that\n"
- "the field is bigger at z = 'Exit Parameter 1' than at\n"
- "z = 'Exit Parameter 3'");
+ if (computeEntranceFringe(entranceParameter1_m) < computeEntranceFringe(entranceParameter3_m)) {
+ for (int index = 0; index < polyOrderEntry_m + 1; ++ index) {
+ if (index % 2 == 0) continue;
+ engeCoeffsEntry_m[index] *= -1;
+ }
+ }
+
+ if (computeExitFringe(exitParameter1_m) < computeExitFringe(exitParameter3_m)) {
+ for (int index = 0; index < polyOrderExit_m + 1; ++ index) {
+ if (index % 2 == 0) continue;
+ engeCoeffsExit_m[index] *= -1;
+ }
+ }
}
void FM1DProfile1::freeMap() {
diff --git a/src/Classic/Utilities/MSLang/Repeat.cpp b/src/Classic/Utilities/MSLang/Repeat.cpp
index 618c8baf5a8a9e2696ed24472fa5611fdf8e3b9d..19781de6b7c4afdf3bbc57bf2aba9c3b86f40a75 100644
--- a/src/Classic/Utilities/MSLang/Repeat.cpp
+++ b/src/Classic/Utilities/MSLang/Repeat.cpp
@@ -2,8 +2,6 @@
#include "Utilities/MSLang/ArgumentExtractor.h"
#include "Utilities/MSLang/matheval.h"
-#include <boost/regex.hpp>
-
namespace mslang {
void Repeat::print(int indentwidth) {
std::string indent(indentwidth, ' ');
@@ -24,7 +22,7 @@ namespace mslang {
const unsigned int size = bfuncs.size();
AffineTransformation current_trafo = trafo;
- for (unsigned int i = 0; i < N_m; ++ i) {
+ for (int i = 0; i < N_m; ++ i) {
for (unsigned int j = 0; j < size; ++ j) {
Base *obj = bfuncs[j]->clone();
obj->trafo_m = obj->trafo_m.mult(current_trafo);
diff --git a/src/Classic/Utilities/MSLang/Repeat.h b/src/Classic/Utilities/MSLang/Repeat.h
index 0b1eab378fe299aac7c99ad49d6b60a0289495c7..2985465382a85f338656bb5af2d0fd7983be3d61 100644
--- a/src/Classic/Utilities/MSLang/Repeat.h
+++ b/src/Classic/Utilities/MSLang/Repeat.h
@@ -6,7 +6,7 @@
namespace mslang {
struct Repeat: public Function {
Function* func_m;
- unsigned int N_m;
+ signed int N_m;
double shiftx_m;
double shifty_m;
double rot_m;
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 0ad95b65df6253378a7f48bf91b886ce07cd3ef0..99e1185f13c2df94f92137c5cc48bb6d5d929c6f 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -1112,8 +1112,8 @@ void Distribution::createDistributionFromFile(size_t numberOfParticles, double m
// Data input file is only read by node 0.
std::ifstream inputFile;
+ std::string fileName = Attributes::getString(itsAttr[Attrib::Distribution::FNAME]);
if (Ippl::myNode() == 0) {
- std::string fileName = Attributes::getString(itsAttr[Attrib::Distribution::FNAME]);
inputFile.open(fileName.c_str());
if (inputFile.fail())
throw OpalException("Distribution::create()",
@@ -1127,36 +1127,38 @@ void Distribution::createDistributionFromFile(size_t numberOfParticles, double m
* We read in the particle information using node zero, but save the particle
* data to each processor in turn.
*/
- int saveProcessor = -1;
+ unsigned int saveProcessor = 0;
pmean_m = 0.0;
size_t numPartsToSend = 0;
unsigned int distributeFrequency = 1000;
- size_t singleDataSize = (sizeof(int) + 6 * sizeof(double));
- size_t dataSize = distributeFrequency * singleDataSize;
+ size_t singleDataSize = (/*sizeof(int) +*/ 6 * sizeof(double));
+ unsigned int dataSize = distributeFrequency * singleDataSize;
std::vector<char> data;
data.reserve(dataSize);
const char* buffer;
- for (size_t particleIndex = 0; particleIndex < numberOfParticlesRead; ++ particleIndex) {
-
- // Read particle.
- Vector_t R(0.0);
- Vector_t P(0.0);
-
- ++ saveProcessor;
- if (saveProcessor >= Ippl::getNodes())
- saveProcessor = 0;
-
- if (Ippl::myNode() == 0) {
- inputFile >> R(0)
- >> P(0)
- >> R(1)
- >> P(1)
- >> R(2)
- >> P(2);
+ if (Ippl::myNode() == 0) {
+ char lineBuffer[1024];
+ unsigned int numParts = 0;
+ while (!inputFile.eof()) {
+ inputFile.getline(lineBuffer, 1024);
+
+ Vector_t R(0.0), P(0.0);
+
+ std::istringstream line(lineBuffer);
+ line >> R(0);
+ if (line.rdstate()) break;
+ line >> P(0);
+ line >> R(1);
+ line >> P(1);
+ line >> R(2);
+ line >> P(2);
+
+ if (saveProcessor >= (unsigned)Ippl::getNodes())
+ saveProcessor = 0;
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
P(0) = converteVToBetaGamma(P(0), massIneV);
@@ -1166,9 +1168,7 @@ void Distribution::createDistributionFromFile(size_t numberOfParticles, double m
pmean_m += P;
- if (saveProcessor > 0) {
- buffer = reinterpret_cast<const char*>(&saveProcessor);
- data.insert(data.end(), buffer, buffer + sizeof(int));
+ if (saveProcessor > 0u) {
buffer = reinterpret_cast<const char*>(&R[0]);
data.insert(data.end(), buffer, buffer + 3 * sizeof(double));
buffer = reinterpret_cast<const char*>(&P[0]);
@@ -1182,66 +1182,67 @@ void Distribution::createDistributionFromFile(size_t numberOfParticles, double m
pyDist_m.push_back(P(1));
pzDist_m.push_back(P(2));
}
- } else {
- if (saveProcessor > 0) {
- ++ numPartsToSend;
- }
- }
- // split distribution of particles onto cores such that each time
- // <distributionFrequency> number of particles are sent
- if (numPartsToSend % distributeFrequency == distributeFrequency - 1) {
- MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
- numPartsToSend = 0;
+ ++ numParts;
+ ++ saveProcessor;
+
+ if (numPartsToSend % distributeFrequency == distributeFrequency - 1) {
+ MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
+ MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
+ numPartsToSend = 0;
- if (Ippl::myNode() == 0) {
std::vector<char>().swap(data);
data.reserve(dataSize);
- } else {
- size_t i = 0;
- while (i < dataSize) {
- int saveProcessor = *reinterpret_cast<const int*>(&data[i]);
-
- if (saveProcessor == Ippl::myNode()) {
- i += sizeof(int);
- const double *tmp = reinterpret_cast<const double*>(&data[i]);
- xDist_m.push_back(tmp[0]);
- yDist_m.push_back(tmp[1]);
- tOrZDist_m.push_back(tmp[2]);
- pxDist_m.push_back(tmp[3]);
- pyDist_m.push_back(tmp[4]);
- pzDist_m.push_back(tmp[5]);
- i += 6 * sizeof(double);
- } else {
- i += singleDataSize;
- }
- }
}
}
- }
- // send remaining particles (less than <distributionFrequency>)
- MPI_Bcast(&data[0], numPartsToSend * singleDataSize, MPI_CHAR, 0, Ippl::getComm());
+ dataSize = (numberOfParticlesRead == numParts? data.size(): std::numeric_limits<unsigned int>::max());
+ MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
+ if (numberOfParticlesRead != numParts) {
+ throw OpalException("Distribution::createDistributionFromFile",
+ "Found " +
+ std::to_string(numParts) +
+ " particles in file '" +
+ fileName +
+ "' instead of " +
+ std::to_string(numberOfParticlesRead));
+ }
+ MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
- if (Ippl::myNode() > 0) {
- size_t i = 0;
- while (i < numPartsToSend * singleDataSize) {
- int saveProcessor = *reinterpret_cast<const int*>(&data[i]);
-
- if (saveProcessor == Ippl::myNode()) {
- i += sizeof(int);
- const double *tmp = reinterpret_cast<const double*>(&data[i]);
- xDist_m.push_back(tmp[0]);
- yDist_m.push_back(tmp[1]);
- tOrZDist_m.push_back(tmp[2]);
- pxDist_m.push_back(tmp[3]);
- pyDist_m.push_back(tmp[4]);
- pzDist_m.push_back(tmp[5]);
- i += 6 * sizeof(double);
- } else {
- i += singleDataSize;
+ } else {
+ do {
+ size_t i = 0;
+ MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
+ if (dataSize == std::numeric_limits<unsigned int>::max()) {
+ throw OpalException("Distribution::createDistributionFromFile",
+ "Couldn't find " +
+ std::to_string(numberOfParticlesRead) +
+ " particles in file '" +
+ fileName + "'");
}
- }
+ MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
+
+ while (i < dataSize) {
+
+ if (saveProcessor + 1 == (unsigned) Ippl::myNode()) {
+ const double *tmp = reinterpret_cast<const double*>(&data[i]);
+ xDist_m.push_back(tmp[0]);
+ yDist_m.push_back(tmp[1]);
+ tOrZDist_m.push_back(tmp[2]);
+ pxDist_m.push_back(tmp[3]);
+ pyDist_m.push_back(tmp[4]);
+ pzDist_m.push_back(tmp[5]);
+ i += 6 * sizeof(double);
+ } else {
+ i += singleDataSize;
+ }
+
+ ++ saveProcessor;
+ if (saveProcessor + 1 >= (unsigned) Ippl::getNodes()) {
+ saveProcessor = 0;
+ }
+ }
+ } while (dataSize == distributeFrequency * singleDataSize);
}
pmean_m /= numberOfParticlesRead;
@@ -1281,20 +1282,20 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
throw OpalException("Distribution::createMatchedGaussDistribution",
"didn't find any Cyclotron element in line");
}
-
+
/* FIXME we need to remove const-ness otherwise we can't update the injection radius
* and injection radial momentum. However, there should be a better solution ..
*/
Cyclotron* CyclotronElement = const_cast<Cyclotron*>(CyclotronVisitor.front());
-
+
bool full = !Attributes::getBool(itsAttr[Attrib::Distribution::SECTOR]);
int Nint = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSTEPS]);
-
+
if ( Nint < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative number of integration steps");
-
+
*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;
@@ -1320,37 +1321,37 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
"Missing attributes 'FMLOWE' and/or 'FMHIHGE' in "
"'CYCLOTRON' definition.");
}
-
-
+
+
std::size_t maxitCOF =
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSCO]);
-
+
double rguess =
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]);
-
+
int nSector = (int)CyclotronElement->getSymmetry();
-
+
double accuracy =
Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]);
-
+
if ( Options::cloTuneOnly ) {
*gmsg << "* Stopping after closed orbit and tune calculation" << endl;
typedef std::vector<double> container_t;
typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
-
+
cof_t cof(E_m*1E-6, massIneV*1E-6, wo, Nint, accuracy,
maxitCOF, fmLowE, fmHighE, nSector,
CyclotronElement->getFieldMapFN(), rguess,
CyclotronElement->getCyclotronType(),
CyclotronElement->getBScale(), false);
-
+
std::pair<double, double> tunes = cof.getTunes();
double ravg = cof.getAverageRadius(); // average radius
-
+
container_t reo = cof.getOrbit(CyclotronElement->getPHIinit());
container_t peo = cof.getMomentum(CyclotronElement->getPHIinit());
-
+
*gmsg << "* ----------------------------" << endl
<< "* Closed orbit info (Gordon units):" << endl
@@ -1362,12 +1363,12 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
<< "* horizontal tune: " << tunes.first << endl
<< "* vertical tune: " << tunes.second << endl
<< "* ----------------------------" << endl << endl;
-
+
std::exit(0);
}
-
+
bool writeMap = true;
-
+
typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
@@ -1385,7 +1386,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
CyclotronElement->getBScale(),
writeMap);
-
+
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
maxitCOF,
@@ -1410,9 +1411,9 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
*gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
- *gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
+ *gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
-
+
}
*gmsg << " \\\\" << endl;
}
@@ -1433,24 +1434,24 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
// change units from mm to m
for (unsigned int i = 0; i < 6; ++ i)
for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
-
+
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.");
+ "Negative value on the diagonal of the sigma matrix.");
}
-
+
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
sigmaR_m[1] = std::sqrt(sigma(4, 4));
-
+
//p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
-
+
double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma -
sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
-
+
double pl = std::sqrt(pl2);
sigmaP_m[1] = gamma*(pl - beta*gamma);
@@ -1463,7 +1464,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
correlationMatrix_m(5, 1) = sigma(1, 5) / (sqrt(sigma(1, 1) * sigma(5, 5)));
createDistributionGauss(numberOfParticles, massIneV);
-
+
// update injection radius and radial momentum
CyclotronElement->setRinit(siggen->getInjectionRadius() * 1.0e3);
CyclotronElement->setPRinit(siggen->getInjectionMomentum());
@@ -2550,7 +2551,7 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
double rn = 1e-12;
while (errcode == GSL_EDOM) {
-
+
// Resets the correlation matrix
for (unsigned int i = 0; i < 6; ++ i) {
gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
@@ -3570,7 +3571,7 @@ void Distribution::setAttributes() {
itsAttr[Attrib::Distribution::NSTEPS]
= Attributes::makeReal("NSTEPS", "Number of integration steps of closed orbit finder (matched gauss)"
" (default: 720)", 720);
-
+
itsAttr[Attrib::Distribution::RGUESS]
= Attributes::makeReal("RGUESS", "Guess value of radius (m) for closed orbit finder ", -1);
@@ -4832,4 +4833,4 @@ void Distribution::adjustPhaseSpace(double massIneV) {
// mode:c++
// c-basic-offset: 4
// indent-tabs-mode:nil
-// End:
+// End:
\ No newline at end of file
diff --git a/src/Track/TrackRun.cpp b/src/Track/TrackRun.cpp
index 45617bc8c7ca0e1fc973374476da17f92dc0884b..f863e52f78752af30720c7566ef6cbd211c8604e 100644
--- a/src/Track/TrackRun.cpp
+++ b/src/Track/TrackRun.cpp
@@ -73,7 +73,6 @@ namespace {
// The attributes of class TrackRun.
enum {
METHOD, // Tracking method to use.
- FNAME, // The name of file to be written.
TURNS, // The number of turns to be tracked.
MBMODE, // The working way for multi-bunch mode for OPAL-cycl: "FORCE" or "AUTO"
PARAMB, // The control parameter for "AUTO" mode of multi-bunch,
@@ -111,14 +110,11 @@ TrackRun::TrackRun():
itsAttr[PARAMB] = Attributes::makeReal
("PARAMB", " Control parameter to define when to start multi-bunch mode, only available in \"AUTO\" mode ", 5.0);
-
+
itsAttr[MB_ETA] = Attributes::makeReal("MB_ETA",
"The scale parameter for binning in multi-bunch mode",
0.01);
- itsAttr[FNAME] = Attributes::makeString
- ("FILE", "Name of file to be written", "TRACK");
-
itsAttr[BEAM] = Attributes::makeString
("BEAM", "Name of beam ", "BEAM");
itsAttr[FIELDSOLVER] = Attributes::makeString
@@ -197,7 +193,7 @@ void TrackRun::execute() {
Track::block->bunch, Track::block->reference,
false, false);
} else if(method == "THICK") {
- setupThickTracker();
+ setupThickTracker();
// } else if(method == "PARALLEL-SLICE" || method == "OPAL-E") {
// setupSliceTracker();
} else if(method == "PARALLEL-T" || method == "OPAL-T") {
@@ -212,15 +208,6 @@ void TrackRun::execute() {
}
if(method == "THIN" || method == "THICK") {
- //
- std::string file = Attributes::getString(itsAttr[FNAME]);
- std::ofstream os(file.c_str());
-
- if(os.bad()) {
- throw OpalException("TrackRun::execute()",
- "Unable to open output file \"" + file + "\".");
- }
-
int turns = int(Round(Attributes::getReal(itsAttr[TURNS])));
// Track for the all but last turn.
@@ -228,7 +215,7 @@ void TrackRun::execute() {
itsTracker->execute();
}
- // Track the last turn.
+ // Track the last turn.
itsTracker->execute();
} else {
@@ -471,8 +458,8 @@ void TrackRun::setupThickTracker()
itsTracker = new ThickTracker(*Track::block->use->fetchLine(),
Track::block->bunch, *ds, Track::block->reference,
- false, false, Track::block->localTimeSteps,
- Track::block->zstart, Track::block->zstop, Track::block->dT,
+ false, false, Track::block->localTimeSteps,
+ Track::block->zstart, Track::block->zstop, Track::block->dT,
Track::block->truncOrder);
}
@@ -824,10 +811,10 @@ void TrackRun::setupCyclotronTracker(){
itsTracker->setMultiBunchMode("FORCE");
}
} else {
- std::string mbmode = Util::toUpper(Attributes::getString(itsAttr[MBMODE]));
+ std::string mbmode = Util::toUpper(Attributes::getString(itsAttr[MBMODE]));
itsTracker->setMultiBunchMode(mbmode);
}
-
+
dynamic_cast<ParallelCyclotronTracker*>(itsTracker)->setMultiBunchEta(Attributes::getReal(itsAttr[MB_ETA]));
}
}