diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 3c7135b14824ef1eeee820cc13895bc48429878b..54f214a1d7253798fc6cdec062d663b30ef8a63e 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -2428,8 +2428,8 @@ void ParallelCyclotronTracker::checkFileMomentum() {
allreduce(pTotalMean, 1, std::plus<double>());
pTotalMean /= initialTotalNum_m;
- double tolerance = itsReference.getMomentumTolerance();
+ double tolerance = itsBunch_m->getMomentumTolerance();
if (tolerance > 0 &&
std::abs(pTotalMean - referencePtot) / pTotalMean > tolerance) {
throw OpalException("ParallelCyclotronTracker::checkFileMomentum",
diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h
index 2ad0bcfb168b0ccb0ad8ccc65a64e26db76e99b5..361286101e80e1378de8b98b28b24c7bc774e5fa 100644
--- a/src/Classic/Algorithms/PartBunchBase.h
+++ b/src/Classic/Algorithms/PartBunchBase.h
@@ -447,6 +447,8 @@ public:
DistributionType getDistType() const;
+ double getMomentumTolerance() const;
+
Quaternion_t getQKs3D();
void setQKs3D(Quaternion_t q);
Vector_t getKs3DRefr();
diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp
index d90e0d3652fca7fcb9bc0bd273153a611823c39c..a4bc0edda636198a13485dee7d6d8ee52bb28230 100644
--- a/src/Classic/Algorithms/PartBunchBase.hpp
+++ b/src/Classic/Algorithms/PartBunchBase.hpp
@@ -1678,6 +1678,12 @@ DistributionType PartBunchBase<T, Dim>::getDistType() const {
}
+template <class T, unsigned Dim>
+double PartBunchBase<T, Dim>::getMomentumTolerance() const {
+ return reference->getMomentumTolerance();
+}
+
+
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::resetQ(double q) {
const_cast<PartData *>(reference)->setQ(q);
diff --git a/src/Classic/Algorithms/PartData.cpp b/src/Classic/Algorithms/PartData.cpp
index a533942157ef29ab73148afb136f0fc898c15666..68c3c9aeec4aac90d8af098e1eb728f72824f883 100644
--- a/src/Classic/Algorithms/PartData.cpp
+++ b/src/Classic/Algorithms/PartData.cpp
@@ -1,95 +1,90 @@
-// ------------------------------------------------------------------------
-// $RCSfile: PartData.cpp,v $
-// ------------------------------------------------------------------------
-// $Revision: 1.1.1.1.2.1 $
-// ------------------------------------------------------------------------
-// Copyright: see Copyright.readme
-// ------------------------------------------------------------------------
//
-// Class: PartData
+// Class PartData
// PartData represents a set of reference values for use in algorithms.
//
-// ------------------------------------------------------------------------
-// Class category: Algorithms
-// ------------------------------------------------------------------------
+// Copyright (c) 200x - 2023, Paul Scherrer Institut, Villigen PSI, Switzerland
+// All rights reserved
//
-// $Date: 2003/12/02 23:04:59 $
-// $Author: dbruhwil $
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-// ------------------------------------------------------------------------
-
#include "Algorithms/PartData.h"
+
#include "Utilities/LogicalError.h"
-#include <cmath>
+#include <cmath>
-// Class PartData
-// ------------------------------------------------------------------------
PartData::PartData(double q, double m, double momentum) {
- charge = q;
- mass = m;
+ charge_m = q;
+ mass_m = m;
setP(momentum);
}
PartData::PartData() {
- charge = 1.0;
- mass = 0.0;
- beta = 1.0;
- gamma = 1.0e10;
+ charge_m = 1.0;
+ mass_m = 0.0;
+ beta_m = 1.0;
+ gamma_m = 1.0e10;
}
void PartData::setP(double p) {
- if(mass == 0.0) {
+ if (mass_m == 0.0) {
throw LogicalError("PartData::setP()",
"Particle mass must not be zero.");
}
- if(p == 0.0) {
+ if (p == 0.0) {
throw LogicalError("PartData::setP()",
"Particle momentum must not be zero.");
}
- double e = std::sqrt(p * p + mass * mass);
- beta = p / e;
- gamma = e / mass;
+ double e = std::sqrt(p * p + mass_m * mass_m);
+ beta_m = p / e;
+ gamma_m = e / mass_m;
}
void PartData::setE(double energy) {
- if(energy <= mass) {
+ if (energy <= mass_m) {
throw LogicalError("PartData::setE()", "Energy should be > mass.");
}
- gamma = energy / mass;
+ gamma_m = energy / mass_m;
//beta = std::sqrt(energy*energy - mass*mass) / energy;
- double ginv = 1.0 / gamma;
- beta = std::sqrt((1.0 - ginv) * (1.0 + ginv));
+ double ginv = 1.0 / gamma_m;
+ beta_m = std::sqrt((1.0 - ginv) * (1.0 + ginv));
}
void PartData::setBeta(double v) {
- if(v >= 1.0) {
+ if (v >= 1.0) {
throw LogicalError("PartData::setBeta()", "Beta should be < 1.");
}
- beta = v;
- gamma = 1.0 / std::sqrt(1.0 - beta * beta);
+ beta_m = v;
+ gamma_m = 1.0 / std::sqrt(1.0 - beta_m * beta_m);
}
void PartData::setGamma(double v) {
- if(v <= 1.0) {
+ if (v <= 1.0) {
throw LogicalError("PartData::setGamma()", "Gamma should be > 1.");
}
- gamma = v;
- beta = std::sqrt(gamma * gamma - 1.0) / gamma;
+ gamma_m = v;
+ beta_m = std::sqrt(gamma_m * gamma_m - 1.0) / gamma_m;
}
-
void PartData::setMomentumTolerance(double tolerance) {
- momentumTolerance = tolerance;
+ momentumTolerance_m = tolerance;
}
diff --git a/src/Classic/Algorithms/PartData.h b/src/Classic/Algorithms/PartData.h
index 90149e29d199cf8b8b1afcf4b983c3d6c78a1df9..d335c1735945a4e7fd7d8c988577b032ed368e26 100644
--- a/src/Classic/Algorithms/PartData.h
+++ b/src/Classic/Algorithms/PartData.h
@@ -1,44 +1,43 @@
-#ifndef MAD_PartData_HH
-#define MAD_PartData_HH
-
-// ------------------------------------------------------------------------
-// $RCSfile: PartData.h,v $
-// ------------------------------------------------------------------------
-// $Revision: 1.1.1.1 $
-// ------------------------------------------------------------------------
-// Copyright: see Copyright.readme
-// ------------------------------------------------------------------------
-// Description:
//
-// ------------------------------------------------------------------------
-// Class category: Algorithms
-// ------------------------------------------------------------------------
+// Class PartData
+// PartData represents a set of reference values for use in algorithms.
//
-// $Date: 2000/03/27 09:32:33 $
-// $Author: fci $
+// Copyright (c) 200x - 2023, Paul Scherrer Institut, Villigen PSI, Switzerland
+// All rights reserved
//
-// ------------------------------------------------------------------------
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
+#ifndef MAD_PartData_HH
+#define MAD_PartData_HH
-// Class PartData
-// ------------------------------------------------------------------------
-/// Particle reference data.
-// This class encapsulates the reference data for a beam:
-// [UL]
-// [LI]charge per particle expressed in proton charges,
-// [LI]mass per particle expressed in eV,
-// [LI]reference momentum per particle expressed in eV.
-// [LI]momentumTolerance Fractional tolerance to deviations in the distribution
-// compared to the reference data at initialisation
-// If negative, no tolerance checking is done.
-// [/UL]
-// The copy constructor, destructor, and assignment operator generated
-// by the compiler perform the correct operation. For speed reasons
-// they are not implemented.
+/** Class PartData
+ * ------------------------------------------------------------------------
+ * Particle reference data.
+ * This class encapsulates the reference data for a beam:
+ * [UL]
+ * [LI]charge per particle expressed in proton charges,
+ * [LI]mass per particle expressed in eV,
+ * [LI]reference momentum per particle expressed in eV.
+ * [LI]momentumTolerance Fractional tolerance to deviations in the distribution
+ * compared to the reference data at initialisation
+ * If negative, no tolerance checking is done.
+ * [/UL]
+ * The copy constructor, destructor, and assignment operator generated
+ * by the compiler perform the correct operation. For speed reasons
+ * they are not implemented.
+ */
class PartData {
public:
-
/// Constructor.
// Inputs are:
// [DL]
@@ -68,6 +67,15 @@ public:
/// The relativistic gamma per particle.
double getGamma() const;
+ /// Get the momentum tolerance
+ double getMomentumTolerance() const;
+
+ /// Set reference charge expressed in proton charges.
+ void setQ(double q);
+
+ /// Set reference mass expressed in eV/c^2.
+ void setM(double m);
+
/// Set reference momentum.
// Input is the momentum in eV.
void setP(double p);
@@ -84,63 +92,55 @@ public:
// Input is the relativistic gamma = E/(m*c*c).
void setGamma(double gamma);
- /// Get the momentum tolerance
- double getMomentumTolerance() const;
-
- /// Set the momentum tolerance
+ /// Set the momentum tolerance.
void setMomentumTolerance(double tolerance);
- /// Set reference mass expressed in eV/c^2
- inline void setM(double m){mass = m;}
-
- /// Set reference charge expressed in proton charges,
- inline void setQ(double q) {charge = q;}
-
protected:
-
// The reference information.
- double charge; // Particle charge.
- double mass; // Particle mass.
- double beta; // particle velocity divided by c.
- double gamma; // particle energy divided by particle mass
- double momentumTolerance = 1e-2; // tolerance to momentum deviations
+ double charge_m; // Particle charge.
+ double mass_m; // Particle mass.
+ double beta_m; // Particle velocity divided by c.
+ double gamma_m; // Particle energy divided by particle mass.
+ double momentumTolerance_m = 1e-2; // Tolerance to momentum deviations.
};
// Inline functions.
// ------------------------------------------------------------------------
+inline void PartData::setM(double m) {
+ mass_m = m;
+}
-inline double PartData::getQ() const {
- return charge;
+inline void PartData::setQ(double q) {
+ charge_m = q;
}
+inline double PartData::getQ() const {
+ return charge_m;
+}
inline double PartData::getM() const {
- return mass;
+ return mass_m;
}
-
inline double PartData::getP() const {
- return beta * gamma * mass;
+ return beta_m * gamma_m * mass_m;
}
-
inline double PartData::getE() const {
- return gamma * mass;
+ return gamma_m * mass_m;
}
-
inline double PartData::getBeta() const {
- return beta;
+ return beta_m;
}
-
inline double PartData::getGamma() const {
- return gamma;
+ return gamma_m;
}
inline double PartData::getMomentumTolerance() const {
- return momentumTolerance;
+ return momentumTolerance_m;
}
#endif // MAD_PartData_HH
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 2fbf3dd24390ee4eae7a82d0f0b9291664d91e96..478f64eb760aaa12d234930f859911df38c31993 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -868,11 +868,11 @@ void Distribution::checkParticleNumber(size_t &numberOfParticles) {
}
}
-void Distribution::checkFileMomentum() {
+void Distribution::checkFileMomentum(double momentumTol) {
// If the distribution was read from a file, the file momentum pmean_m[2]
// should coincide with the momentum given in the beam command avrgpz_m.
-
- if (std::abs(pmean_m[2] - avrgpz_m) / pmean_m[2] > 1e-2) {
+ if (momentumTol > 0 &&
+ std::abs(pmean_m[2] - avrgpz_m) / pmean_m[2] > momentumTol) {
throw OpalException("Distribution::checkFileMomentum",
"The z-momentum of the particle distribution\n" +
std::to_string(pmean_m[2]) + "\n"
@@ -1479,7 +1479,8 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
checkEmissionParameters();
} else {
if (distrTypeT_m == DistributionType::FROMFILE) {
- checkFileMomentum();
+ double momentumTol = beam->getMomentumTolerance();
+ checkFileMomentum(momentumTol);
} else {
pmean_m = Vector_t(0, 0, avrgpz_m);
}
@@ -2925,7 +2926,6 @@ void Distribution::printDist(Inform &os, size_t numberOfParticles) const {
throw OpalException("Distribution::printDist",
"Unknown \"TYPE\" of \"DISTRIBUTION\"");
}
-
}
void Distribution::printDistBinomial(Inform &os) const {
@@ -3337,8 +3337,7 @@ void Distribution::setAttributes() {
itsAttr[Attrib::Distribution::FNAME]
= Attributes::makeString("FNAME", "File for reading in 6D particle "
- "coordinates.", "");
-
+ "coordinates.");
itsAttr[Attrib::Distribution::WRITETOFILE]
= Attributes::makeBool("WRITETOFILE", "Write initial distribution to file.",
diff --git a/src/Distribution/Distribution.h b/src/Distribution/Distribution.h
index 643df675bde14ca6b213eefa20c3122f6009e1c9..4c9e66c5df975ed01027116d5f364b095859afe8 100644
--- a/src/Distribution/Distribution.h
+++ b/src/Distribution/Distribution.h
@@ -295,7 +295,7 @@ private:
void checkEmissionParameters();
void checkIfEmitted();
void checkParticleNumber(size_t &numberOfParticles);
- void checkFileMomentum();
+ void checkFileMomentum(double momentumTol);
void chooseInputMomentumUnits(InputMomentumUnits inputMoUnits);
size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
diff --git a/src/PyOpal/PyObjects/PyBeam.cpp b/src/PyOpal/PyObjects/PyBeam.cpp
index 73b1c0bc08ea67456cd5c2423a2f54d2c25f6299..5f671159f15b681e5105fdba108a0fd3b32e889e 100644
--- a/src/PyOpal/PyObjects/PyBeam.cpp
+++ b/src/PyOpal/PyObjects/PyBeam.cpp
@@ -35,7 +35,7 @@ std::vector<PyOpalObjectNS::AttributeDef> PyOpalObjectNS::PyOpalObject<Beam>::at
{"BCURRENT", "beam_current", "", PyOpalObjectNS::DOUBLE},
{"BFREQ", "beam_frequency", "", PyOpalObjectNS::DOUBLE},
{"NPART", "number_of_particles", "", PyOpalObjectNS::DOUBLE},
- {"MOMENTUM_TOLERANCE", "momentum_tolerance", "", PyOpalObjectNS::DOUBLE},
+ {"MOMENTUMTOLERANCE", "momentum_tolerance", "", PyOpalObjectNS::DOUBLE},
};
BOOST_PYTHON_MODULE(beam) {
diff --git a/src/Structure/Beam.cpp b/src/Structure/Beam.cpp
index 343f1b22141823e65193f2613754b9c97437247d..c3dcc1115c4fbf4cbaff4f02ad59cd1b4bc752d3 100644
--- a/src/Structure/Beam.cpp
+++ b/src/Structure/Beam.cpp
@@ -4,7 +4,7 @@
// A BEAM definition is used by most physics commands to define the
// particle charge and the reference momentum, together with some other data.
//
-// Copyright (c) 200x - 2021, Paul Scherrer Institut, Villigen PSI, Switzerland
+// Copyright (c) 200x - 2023, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
@@ -46,7 +46,7 @@ namespace {
BCURRENT, // Beam current in A
BFREQ, // Beam frequency in MHz
NPART, // Number of particles per bunch
- MOMENTUM_TOLERANCE, // Fractional tolerance to momentum deviations
+ MOMENTUMTOLERANCE, // Fractional tolerance to momentum deviations
SIZE
};
}
@@ -97,11 +97,11 @@ Beam::Beam():
itsAttr[NPART] = Attributes::makeReal
("NPART", "Number of particles in bunch");
- itsAttr[MOMENTUM_TOLERANCE] = Attributes::makeReal
- ("MOMENTUM_TOLERANCE",
- "Fractional tolerance to deviations in the distribution compared to the"
- "reference data at initialisation. If <= 0, no tolerance checking is"
- "done (default=1e-2).", 1e-2);
+ itsAttr[MOMENTUMTOLERANCE] = Attributes::makeReal
+ ("MOMENTUMTOLERANCE",
+ "Fractional tolerance to deviations in the distribution "
+ "compared to the reference data at initialisation (default=1e-2). "
+ "If MOMENTUMTOLERANCE<=0, no tolerance checking is done.", 1e-2);
// Set up default beam.
Beam* defBeam = clone("UNNAMED_BEAM");
@@ -205,10 +205,6 @@ double Beam::getFrequency() const {
return Attributes::getReal(itsAttr[BFREQ]);
}
-double Beam::getMomentumTolerance() const {
- return Attributes::getReal(itsAttr[MOMENTUM_TOLERANCE]);
-}
-
double Beam::getChargePerParticle() const {
return std::copysign(1.0, getCharge()) * getCurrent()
/ (getFrequency() * Units::MHz2Hz)
@@ -259,8 +255,10 @@ void Beam::update() {
"\"PC\" should be greater than 0.");
}
}
- double tol = Attributes::getReal(itsAttr[MOMENTUM_TOLERANCE]);
- reference.setMomentumTolerance(tol);
+
+ double momentumTol = Attributes::getReal(itsAttr[MOMENTUMTOLERANCE]);
+ reference.setMomentumTolerance(momentumTol);
+
// Set default name.
if (getOpalName().empty()) setOpalName("UNNAMED_BEAM");
}
@@ -273,8 +271,11 @@ void Beam::print(std::ostream& os) const {
<< "* PARTICLE " << Attributes::getString(itsAttr[PARTICLE]) << '\n'
<< "* REST MASS " << Attributes::getReal(itsAttr[MASS]) << " [GeV]\n"
<< "* CHARGE " << (charge > 0 ? '+' : '-') << "e * " << std::abs(charge) << " \n"
- << "* MOMENTUM " << reference.getP() << " [eV/c]\n"
- << "* CURRENT " << Attributes::getReal(itsAttr[BCURRENT]) << " [A]\n"
+ << "* MOMENTUM " << reference.getP() << " [eV/c]\n";
+ if (Attributes::getReal(itsAttr[MOMENTUMTOLERANCE]) > 0) {
+ os << "* MOMENTUM TOLERANCE " << Attributes::getReal(itsAttr[MOMENTUMTOLERANCE]) << '\n';
+ }
+ os << "* CURRENT " << Attributes::getReal(itsAttr[BCURRENT]) << " [A]\n"
<< "* FREQUENCY " << Attributes::getReal(itsAttr[BFREQ]) << " [MHz]\n"
<< "* NPART " << Attributes::getReal(itsAttr[NPART]) << '\n';
os << "* ********************************************************************************** " << std::endl;
diff --git a/src/Structure/Beam.h b/src/Structure/Beam.h
index 00dc71126742bcc5d728bbb5757ae4c870d28c9a..4e2640b5c8df94686d6014ff74f5d11113ccc44f 100644
--- a/src/Structure/Beam.h
+++ b/src/Structure/Beam.h
@@ -4,7 +4,7 @@
// A BEAM definition is used by most physics commands to define the
// particle charge and the reference momentum, together with some other data.
//
-// Copyright (c) 200x - 2021, Paul Scherrer Institut, Villigen PSI, Switzerland
+// Copyright (c) 200x - 2023, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
@@ -70,9 +70,6 @@ public:
/// Return Particle's rest mass in GeV
double getMass() const;
- /// Return the fractional momentum tolerance
- double getMomentumTolerance() const;
-
/// Charge per macro particle in C
double getChargePerParticle() const;
@@ -94,9 +91,6 @@ private:
// The particle reference data.
PartData reference;
-
- // The converstion from GeV to eV.
- static const double energy_scale;
};
inline std::ostream &operator<<(std::ostream& os, const Beam& b) {