19#include "BDSBunchHaloFlatSigma.hh"
21#include "BDSException.hh"
22#include "BDSParticleCoordsFull.hh"
24#include "parser/beam.h"
28#include "CLHEP/Random/RandFlat.h"
29#include "CLHEP/Units/SystemOfUnits.h"
30#include "Randomize.hh"
43 G4double action = aa.action;
44 G4double angle = aa.angle;
45 G4double alpha = tp.alpha;
46 G4double beta = tp.beta;
47 G4double x = std::sqrt(2 * action * beta) * std::cos(angle);
48 G4double xp = (-std::sqrt(2 * action / beta) * (std::sin(angle) + alpha * std::cos(angle)));
53 EllipsePointGenerator::EllipsePointGenerator(G4double actionIn,
55 CLHEP::RandFlat* flatRandomGeneratorIn):
58 flatRandomGenerator(flatRandomGeneratorIn)
63 std::vector<PhaseSpaceCoord> points;
64 points.reserve(npoints);
67 for (G4int i = 0; i <= npoints; ++i)
69 G4double angle = i * CLHEP::twopi / npoints;
70 angles.push_back(angle);
72 points.push_back(PhaseSpaceCoordFromActionAngle(aa, twisspair));
78 std::vector<G4double> cumulativeDistances = CumulativeDistances(points);
79 cumulativeDistances.insert(std::begin(cumulativeDistances), 0);
80 pathLengths = std::move(cumulativeDistances);
86 G4double pathLength = EllipsePerimeter() * flatRandomGenerator->shoot();
88 G4double angle = PathLengthToAngle(pathLength);
90 return PhaseSpaceCoordFromActionAngle(aa, twisspair);
93 G4double EllipsePointGenerator::PathLengthToAngle(G4double pathLength)
const
96 auto it = std::lower_bound(std::begin(pathLengths), std::end(pathLengths), pathLength);
97 auto n = std::distance(std::begin(pathLengths), it);
102 G4double s0 = *(--it);
103 G4double angle0 = *std::next(std::begin(angles), n - 1);
104 G4double angle1 = *std::next(std::begin(angles), n);
107 G4double angle = angle0 + (pathLength - s0) * ((angle1 - angle0) / (s1 - s0));
112BDSBunchHaloFlatSigma::BDSBunchHaloFlatSigma():
114 alphaX(0.0), alphaY(0.0),
115 betaX(0.0), betaY(0.0),
116 emitX(0.0), emitY(0.0),
117 sigmaX(0.0), sigmaY(0.0),
118 haloNSigmaXInner(0.0), haloNSigmaXOuter(0.0),
119 haloNSigmaYInner(0.0), haloNSigmaYOuter(0.0)
121 flatGen =
new CLHEP::RandFlat(*CLHEP::HepRandom::getTheEngine());
127 G4Transform3D beamlineTransformIn,
128 const G4double beamlineSIn)
152 G4double xsigrange = (haloNSigmaXOuter - haloNSigmaXInner);
153 G4double ysigrange = (haloNSigmaYOuter - haloNSigmaYInner);
154 G4double xnsig = haloNSigmaXInner + flatGen->shoot() * xsigrange;
155 G4double ynsig = haloNSigmaYInner + flatGen->shoot() * ysigrange;
159 G4double csix = std::pow(xnsig *
sigmaX, 2) /
betaX;
160 G4double csiy = std::pow(ynsig *
sigmaY, 2) /
betaY;
162 G4double jx = csix / 2.0;
163 G4double jy = csiy / 2.0;
170 auto xps = epgx.GetRandomPointOnEllipse();
171 auto yps = epgy.GetRandomPointOnEllipse();
173 G4double x = (xps.position +
X0) * CLHEP::m;
174 G4double xp = (xps.momentum +
Xp0) * CLHEP::rad;
175 G4double y = (yps.position +
Y0) * CLHEP::m;
176 G4double yp = (yps.momentum +
Yp0) * CLHEP::rad;
179 G4double weight = 1.0;
189 {
throw BDSException(__METHOD_NAME__,
"emitx must be > 0!");}
191 {
throw BDSException(__METHOD_NAME__,
"emity must be > 0!");}
193 if (haloNSigmaXInner <= 0)
194 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner <= 0");}
196 if (haloNSigmaYInner <= 0)
197 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner <= 0");}
199 if (haloNSigmaXInner > haloNSigmaXOuter)
200 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner cannot be less than haloNSigmaXOuter");}
202 if (haloNSigmaYInner > haloNSigmaYOuter)
203 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner cannot be less than haloNSigmaYOuter");}
G4double alphaX
Twiss parameter.
G4double sigmaY
Twiss parameter.
G4double betaY
Twiss parameter.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
G4double sigmaX
Twiss parameter.
G4double emitX
Twiss parameter.
G4double betaX
Twiss parameter.
G4double emitY
Twiss parameter.
virtual BDSParticleCoordsFull GetNextParticleLocal()
G4double alphaY
Twiss parameter.
virtual void CheckParameters()
The base class for bunch distribution generators.
G4double Yp0
Centre of distributions.
G4double T0
Centre of distributions.
G4double S0
Centre of distributions.
G4double X0
Centre of distributions.
G4double Zp0
Centre of distributions.
G4double Xp0
Centre of distributions.
G4double E0
Centre of distributions.
static void SetEmittances(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, G4double &emittGeometricX, G4double &emittGeometricY, G4double &emittNormalisedX, G4double &emittNormalisedY)
G4double Y0
Centre of distributions.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
static G4double CalculateZp(G4double xp, G4double yp, G4double Zp0)
Calculate zp safely based on other components.
virtual void CheckParameters()
General exception with possible name of object and message.
A set of particle coordinates including energy and weight.
Wrapper for particle definition.
Improve type-safety of native enum data type in C++.
Class for generating points uniformly on ellipse perimeters via interpolation. Part of implementation...
double alfy
initial twiss parameters
double haloNSigmaYInner
for the halo distribution
double betx
initial twiss parameters
double haloNSigmaXInner
for the halo distribution
double haloNSigmaYOuter
for the halo distribution
double haloNSigmaXOuter
for the halo distribution
double alfx
initial twiss parameters
double bety
initial twiss parameters
Return either G4Tubs or G4CutTubs depending on flat face.
Simple struct for storing action/angle pairs to aid readability. Implementation detail.
Simple struct for storing position/momentum pairs to aid readability. Implementation detail.
Simple struct for storing Twiss alpha/beta pairs to aid readability. Implementation detail.