diff --git a/src/Distribution/ClosedOrbitFinder.h b/src/Distribution/ClosedOrbitFinder.h
index 7a351afe91e4690d26ddf87c010c05c72651d8b9..aede80ce4b02f9c84bbb879026a199baee14ec69 100644
--- a/src/Distribution/ClosedOrbitFinder.h
+++ b/src/Distribution/ClosedOrbitFinder.h
@@ -70,6 +70,7 @@ class ClosedOrbitFinder
/// Sets the initial values for the integration and calls findOrbit().
/*!
* @param E0 is the potential energy (particle energy at rest) [MeV].
+ * @param q is the particle charge [e]
* @param N specifies the number of splits (2pi/N), i.e number of integration steps
* @param cycl is the cyclotron element
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
@@ -77,7 +78,7 @@ class ClosedOrbitFinder
* @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over
* in order to avoid first harmonic. Only valid if domain is false
*/
- ClosedOrbitFinder(value_type E0, size_type N,
+ ClosedOrbitFinder(value_type E0, value_type q, size_type N,
Cyclotron* cycl, bool domain = true, int Nsectors = 1);
/// Returns the inverse bending radius (size of container N+1)
@@ -154,7 +155,7 @@ class ClosedOrbitFinder
// void computeVerticalOscillations();
/// This function rotates the calculated closed orbit finder properties to the initial angle
- container_type rotate(value_type angle, container_type& orbitProperty);
+ container_type rotate(value_type angle, const container_type& orbitProperty);
/// Stores current position in horizontal direction for the solutions of the ODE with different initial values
std::array<value_type,2> x_m; // x_m = [x1, x2]
@@ -188,6 +189,9 @@ class ClosedOrbitFinder
/// Is the rest mass [MeV / c**2]
value_type E0_m;
+ // Is the particle charge [e]
+ value_type q_m;
+
/// Is the nominal orbital frequency
/* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
* Coupling by Space Charge in Cyclotrons" (2012), formula (1))
@@ -243,8 +247,8 @@ class ClosedOrbitFinder
/*!
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
*/
- std::function<double(double, double)> bcon_m = [](double e0, double wo) {
- return e0 * 1.0e7 / (/* physics::q0 */ 1.0 * Physics::c * Physics::c / wo);
+ std::function<double(double, double)> bcon_m = [this](double e0, double wo) {
+ return e0 * 1.0e7 / (q_m * Physics::c * Physics::c / wo);
};
Cyclotron* cycl_m;
@@ -258,11 +262,13 @@ template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type,
Size_type,
Stepper>::ClosedOrbitFinder(value_type E0,
+ value_type q,
size_type N, Cyclotron* cycl,
bool domain, int nSectors)
: nxcross_m(0)
, nzcross_m(0)
, E0_m(E0)
+ , q_m(q)
, wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, N_m(N)
, dtheta_m(Physics::two_pi/value_type(N))
@@ -894,15 +900,21 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties(c
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
-ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, container_type &orbitProperty) {
+ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, const container_type &orbitProperty) {
container_type orbitPropertyCopy = orbitProperty;
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
+ if (angle < 0.0) {
+ angle = 360.0 + angle;
+ }
+
// compute starting point
- size_type start = deg_step * angle;
+ unsigned int start = deg_step * angle;
+
+ start %= orbitProperty.size();
// copy end to start
std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 53f20fb0116e25c6a0200b3257c5a5b8622dbc18..1cf37f4dbaa88696e7dc0fab1833f8c3b3679aa7 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -248,7 +248,7 @@ void Distribution::execute() {
void Distribution::update() {
}
-void Distribution::create(size_t &numberOfParticles, double massIneV) {
+void Distribution::create(size_t &numberOfParticles, double massIneV, double charge) {
/*
* If Options::cZero is true than we reflect generated distribution
@@ -284,7 +284,7 @@ void Distribution::create(size_t &numberOfParticles, double massIneV) {
switch (distrTypeT_m) {
case DistrTypeT::MATCHEDGAUSS:
- createMatchedGaussDistribution(numberOfLocalParticles, massIneV);
+ createMatchedGaussDistribution(numberOfLocalParticles, massIneV, charge);
break;
case DistrTypeT::FROMFILE:
createDistributionFromFile(numberOfParticles, massIneV);
@@ -1152,7 +1152,10 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
}
-void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, double massIneV) {
+void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
+ double massIneV,
+ double charge)
+{
/*
ToDo:
@@ -1250,7 +1253,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
- cof_t cof(massIneV*1E-6, Nint, CyclotronElement, full, Nsectors);
+ cof_t cof(massIneV*1E-6, charge, Nint, CyclotronElement, full, Nsectors);
cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
throw EarlyLeaveException("Distribution::CreateMatchedGaussDistribution()",
@@ -1266,6 +1269,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
+ charge,
CyclotronElement,
Nint,
Nsectors,
@@ -1373,7 +1377,7 @@ void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam,
emitting_m = false;
// Create distribution.
- create(numberOfPartToCreate, beam->getM());
+ create(numberOfPartToCreate, beam->getM(), beam->getQ());
// this variable is needed to be compatible with OPAL-T
particlesPerDist_m.push_back(tOrZDist_m.size());
@@ -1453,11 +1457,11 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
setupEmissionModel(beam);
// Create distributions.
- create(particlesPerDist_m.at(0), beam->getM());
+ create(particlesPerDist_m.at(0), beam->getM(), beam->getQ());
size_t distCount = 1;
for (Distribution* addedDist : addedDistributions_m) {
- addedDist->create(particlesPerDist_m.at(distCount), beam->getM());
+ addedDist->create(particlesPerDist_m.at(distCount), beam->getM(), beam->getQ());
distCount++;
}
diff --git a/src/Distribution/Distribution.h b/src/Distribution/Distribution.h
index 31817e3c968886ab2d141f20d0e3c6876441a383..3f963555bc9578ba1fa956bc460a8e3a4d6a9702 100644
--- a/src/Distribution/Distribution.h
+++ b/src/Distribution/Distribution.h
@@ -286,7 +286,13 @@ private:
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelNone(double &pz);
void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
- void create(size_t &numberOfParticles, double massIneV);
+ /*!
+ * Create the particle distribution.
+ * @param numberOfParticles to create
+ * @param massIneV particle charge in eV
+ * @param charge of the particle type in elementary charge
+ */
+ void create(size_t &numberOfParticles, double massIneV, double charge);
void calcPartPerDist(size_t numberOfParticles);
void checkEmissionParameters();
void checkIfEmitted();
@@ -326,7 +332,8 @@ private:
void createDistributionMultiGauss(size_t numberOfParticles, double massIneV);
void createDistributionFromFile(size_t numberOfParticles, double massIneV);
void createDistributionGauss(size_t numberOfParticles, double massIneV);
- void createMatchedGaussDistribution(size_t numberOfParticles, double massIneV);
+ void createMatchedGaussDistribution(size_t numberOfParticles,
+ double massIneV, double charge);
void sampleUniformDisk(gsl_qrng* quasiRandGen2D, double& x1, double& x2);
void fillEBinHistogram();
void fillParticleBins();
diff --git a/src/Distribution/SigmaGenerator.cpp b/src/Distribution/SigmaGenerator.cpp
index 57ee37f2f51467eaf0c50214d792d26095f076fa..cc70c75f17b376807f20d444c784b5cd224bb1df 100644
--- a/src/Distribution/SigmaGenerator.cpp
+++ b/src/Distribution/SigmaGenerator.cpp
@@ -69,6 +69,7 @@ SigmaGenerator::SigmaGenerator(double I,
double ez,
double E,
double m,
+ double q,
const Cyclotron* cycl,
unsigned int N,
unsigned int Nsectors,
@@ -82,6 +83,7 @@ SigmaGenerator::SigmaGenerator(double I,
, nh_m(cycl->getCyclHarm())
, beta_m(std::sqrt(1.0-1.0/gamma2_m))
, m_m(m)
+ , q_m(q)
, niterations_m(0)
, converged_m(false)
, Emin_m(cycl->getFMLowE())
@@ -131,7 +133,7 @@ SigmaGenerator::SigmaGenerator(double I,
// formula (57)
double lam = Physics::two_pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
// [K3] = m
- double K3 = 3.0 * I_m * lam
+ double K3 = 3.0 * std::abs(q_m) * I_m * lam
/ (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m
* Physics::c * beta_m * beta_m * gamma_m * gamma2_m);
@@ -193,7 +195,7 @@ bool SigmaGenerator::match(double accuracy,
std::vector<matrix_t> Mcycs(nSteps_m), Mscs(nSteps_m);
ClosedOrbitFinder<double, unsigned int,
- boost::numeric::odeint::runge_kutta4<container_t> > cof(m_m, N_m, cycl, false, nSectors_m);
+ boost::numeric::odeint::runge_kutta4<container_t> > cof(m_m, q_m, N_m, cycl, false, nSectors_m);
if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
throw OpalException("SigmaGenerator::match()",
@@ -519,13 +521,13 @@ void SigmaGenerator::initialize(double nuz, double ravg)
double lam = Physics::two_pi * Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
// m * c^3 --> c^2 in [m] = eV / c^2 cancel --> m * c in denominator
- double K3 = 3.0 * I_m * lam
+ double K3 = 3.0 * std::abs(q_m) * I_m * lam
/ (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m
* Physics::c * beta_m * beta_m * gamma2_m * gamma_m); // [K3] = m
// c in denominator cancels with unit of [m] = eV / c^2 --> we need to multiply
// with c in order to get dimensionless quantity
- double alpha = Physics::mu_0 * I_m * Physics::c / (5.0 * std::sqrt(10.0) * m
+ double alpha = std::abs(q_m) * Physics::mu_0 * I_m * Physics::c / (5.0 * std::sqrt(10.0) * m
* gamma_m * nh_m) * std::sqrt(rcyc * rcyc * rcyc / (e * e * e)); // [alpha] = 1/rad --> [alpha] = 1
double sig0 = std::sqrt(2.0 * rcyc * e) / gamma_m; // [sig0] = m*sqrt(rad) --> [sig0] = m
diff --git a/src/Distribution/SigmaGenerator.h b/src/Distribution/SigmaGenerator.h
index dc2035047dcf766b3d6c010534ae5311408952c3..f70f034dc7d16c9b2a3793b810bac4a2384547d4 100644
--- a/src/Distribution/SigmaGenerator.h
+++ b/src/Distribution/SigmaGenerator.h
@@ -71,6 +71,7 @@ public:
* @param ez is the emittance in z-direction (vertical), \f$ \left[\varepsilon_{z}\right] = \pi\ mm\ mrad \f$
* @param E is the energy, \f$ \left[E\right] = MeV \f$
* @param m is the mass of the particles \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
+ * @param q is the particle charge [e]
* @param cycl is the cyclotron element
* @param N is the number of integration steps (closed orbit computation). That's why its also the number
* of maps (for each integration step a map)
@@ -79,7 +80,7 @@ public:
* @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
*/
SigmaGenerator(double I, double ex, double ey, double ez,
- double E, double m, const Cyclotron* cycl,
+ double E, double m, double q, const Cyclotron* cycl,
unsigned int N, unsigned int Nsectors, unsigned int truncOrder, bool write = true);
/// Searches for a matched distribution.
@@ -160,6 +161,8 @@ private:
double beta_m;
/// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
double m_m;
+ /// Is the particle charge [e]
+ double q_m;
/// Is the number of iterations needed for convergence
unsigned int niterations_m;
/// Is true if converged, false otherwise