19#include "BDSBunchHalo.hh"
21#include "BDSException.hh"
22#include "BDSParticleCoordsFull.hh"
24#include "parser/beam.h"
28#include "Randomize.hh"
29#include "CLHEP/Random/RandFlat.h"
30#include "CLHEP/Units/SystemOfUnits.h"
34BDSBunchHalo::BDSBunchHalo():
36 alphaX(0.0), alphaY(0.0),
37 betaX(0.0), betaY(0.0),
38 emitX(0.0), emitY(0.0),
39 gammaX(0.0), gammaY(0.0),
40 sigmaX(0.0), sigmaY(0.0),
41 sigmaXp(0.0), sigmaYp(0.0),
42 haloNSigmaXInner(0.0), haloNSigmaXOuter(0.0),
43 haloNSigmaYInner(0.0), haloNSigmaYOuter(0.0),
52 haloPSWeightParameter(0.0),
54 emitInnerX(0.0), emitInnerY(0.0),
55 emitOuterX(0.0), emitOuterY(0.0),
57 xpMax(0.0), ypMax(0.0)
59 flatGen =
new CLHEP::RandFlat(*CLHEP::HepRandom::getTheEngine());
62BDSBunchHalo::~BDSBunchHalo()
70 G4Transform3D beamlineTransformIn,
71 const G4double beamlineSIn)
103 emitInnerX = std::pow(haloNSigmaXInner, 2) *
emitX;
104 emitInnerY = std::pow(haloNSigmaYInner, 2) *
emitY;
105 emitOuterX = std::pow(haloNSigmaXOuter, 2) *
emitX;
106 emitOuterY = std::pow(haloNSigmaYOuter, 2) *
emitY;
108 xMax = haloNSigmaXOuter *
sigmaX;
109 yMax = haloNSigmaYOuter *
sigmaY;
110 xpMax = std::sqrt(std::pow(haloNSigmaXOuter, 2) *
emitX *
gammaX);
111 ypMax = std::sqrt(std::pow(haloNSigmaYOuter, 2) *
emitY *
gammaY);
130 G4double dx = xMax * (1 - 2 * flatGen->shoot());
131 G4double dy = yMax * (1 - 2 * flatGen->shoot());
132 G4double dxp = xpMax * (1 - 2 * flatGen->shoot());
133 G4double dyp = ypMax * (1 - 2 * flatGen->shoot());
136 double emitXSp =
gammaX * std::pow(std::abs(dx), 2) + (2. *
alphaX * dx * dxp) +
betaX * std::pow(std::abs(dxp), 2);
137 double emitYSp =
gammaY * std::pow(std::abs(dy), 2) + (2. *
alphaY * dy * dyp) +
betaY * std::pow(std::abs(dyp), 2);
141 if ((std::abs(emitXSp) < emitInnerX || std::abs(emitYSp) < emitInnerY) ||
142 (std::abs(emitXSp) > emitOuterX || std::abs(emitYSp) > emitOuterY) ||
143 (std::abs(dx) < (haloXCutInner *
sigmaX)) ||
144 (std::abs(dy) < (haloYCutInner *
sigmaY)) ||
145 (std::abs(dx) > (haloXCutOuter *
sigmaX)) ||
146 (std::abs(dy) > (haloYCutOuter *
sigmaY)) ||
147 (std::abs(dxp) < (haloXpCutInner *
sigmaXp)) ||
148 (std::abs(dyp) < (haloYpCutInner *
sigmaYp)) ||
149 (std::abs(dxp) > (haloXpCutOuter *
sigmaXp)) ||
150 (std::abs(dyp) > (haloYpCutOuter *
sigmaYp)) )
159 if (weightFunction ==
"flat" || weightFunction ==
"" || weightFunction ==
"one")
164 else if (weightFunction ==
"oneoverr")
167 wx = std::pow(std::abs(emitInnerX / emitXSp), haloPSWeightParameter);
168 wy = std::pow(std::abs(emitInnerY / emitYSp), haloPSWeightParameter);
170 else if (weightFunction ==
"oneoverrsqrd")
173 double eXsqrd = std::pow(std::abs(emitXSp), 2);
174 double eYsqrd = std::pow(std::abs(emitYSp), 2);
175 double eXInsq = std::pow(std::abs(emitInnerX), 2);
176 double eYInsq = std::pow(std::abs(emitInnerY), 2);
177 wx = std::pow(std::abs(eXInsq / eXsqrd), haloPSWeightParameter);
178 wy = std::pow(std::abs(eYInsq / eYsqrd), haloPSWeightParameter);
180 else if (weightFunction ==
"exp")
182 wx = std::exp(-(emitXSp * haloPSWeightParameter) / (emitInnerX));
183 wy = std::exp(-(emitYSp * haloPSWeightParameter) / (emitInnerY));
187 G4cout << __METHOD_NAME__ << emitXSp/
emitX <<
" " << emitYSp/
emitY <<
" " << wx <<
" " << wy << G4endl;
190 if (flatGen->shoot() > wx && flatGen->shoot() > wy)
196 xp += dxp * CLHEP::rad;
197 yp += dyp * CLHEP::rad;
202 G4cout << __METHOD_NAME__ <<
"selected> " << dx <<
" " << dy <<
" " << dxp <<
" " << dyp << G4endl;
216 {
throw BDSException(__METHOD_NAME__,
"emitx must be finite!");}
218 {
throw BDSException(__METHOD_NAME__,
"emity must be finite!");}
221 {
throw BDSException(__METHOD_NAME__,
"betx must be finite!");}
223 {
throw BDSException(__METHOD_NAME__,
"bety must be finite!");}
225 std::vector<G4String> weightFunctions = {
"",
"one",
"flat",
"oneoverr",
"oneoverrsqrd",
"exp"};
226 auto search = std::find(weightFunctions.begin(), weightFunctions.end(), weightFunction);
227 if (search == weightFunctions.end())
229 G4cerr << __METHOD_END__ <<
"invalid haloPSWeightFunction \"" << weightFunction <<
"\"" << G4endl;
230 G4cout <<
"Available weight functions are:" << G4endl;
231 for (
const auto& w : weightFunctions)
232 {G4cout << w << G4endl;}
236 if (haloNSigmaXInner <= 0)
237 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner <= 0");}
239 if (haloNSigmaYInner <= 0)
240 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner <= 0");}
242 if (haloNSigmaXInner > haloNSigmaXOuter)
243 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner cannot be less than haloNSigmaXOuter");}
245 if (haloNSigmaYInner > haloNSigmaYOuter)
246 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner cannot be less than haloNSigmaYOuter");}
248 if (haloXCutInner < 0)
249 {
throw BDSException(__METHOD_NAME__,
"haloXCutInner < 0");}
251 if (haloYCutInner < 0)
252 {
throw BDSException(__METHOD_NAME__,
"haloYCutInner < 0");}
254 if (haloXCutOuter <= haloXCutInner)
255 {
throw BDSException(__METHOD_NAME__,
"haloXCutOuter must be greater than haloXCutInner!");}
257 if (haloYCutOuter <= haloYCutInner)
258 {
throw BDSException(__METHOD_NAME__,
"haloYCutOuter must be greater than haloYCutInner!");}
260 if (haloXpCutInner < 0)
261 {
throw BDSException(__METHOD_NAME__,
"haloXpCutInner < 0");}
263 if (haloYpCutInner < 0)
264 {
throw BDSException(__METHOD_NAME__,
"haloYpCutInner < 0");}
266 if (haloXpCutOuter <= haloXpCutInner)
267 {
throw BDSException(__METHOD_NAME__,
"haloXpCutOuter must be greater than haloXpCutInner!");}
269 if (haloYpCutOuter <= haloYpCutInner)
270 {
throw BDSException(__METHOD_NAME__,
"haloYpCutOuter must be greater than haloYpCutInner!");}
G4double sigmaX
Twiss parameter.
G4double sigmaYp
Twiss parameter.
G4double alphaY
Twiss parameter.
G4double gammaY
Twiss parameter.
G4double emitX
Twiss parameter.
G4double gammaX
Twiss parameter.
G4double sigmaXp
Twiss parameter.
G4double emitY
Twiss parameter.
G4double sigmaY
Twiss parameter.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
virtual BDSParticleCoordsFull GetNextParticleLocal()
G4double betaY
Twiss parameter.
virtual void CheckParameters()
G4double alphaX
Twiss parameter.
G4double betaX
Twiss parameter.
The base class for bunch distribution generators.
G4double Yp0
Centre of distributions.
G4double T0
Centre of distributions.
G4double S0
Centre of distributions.
G4double Z0
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++.
double alfy
initial twiss parameters
double haloXCutOuter
for the halo distribution
double haloNSigmaYInner
for the halo distribution
double haloXpCutOuter
for the halo distribution
double haloXCutInner
for the halo distribution
double haloYpCutInner
for the halo distribution
double betx
initial twiss parameters
double haloYCutOuter
for the halo distribution
double haloNSigmaXInner
for the halo distribution
double haloNSigmaYOuter
for the halo distribution
double haloYCutInner
for the halo distribution
double haloNSigmaXOuter
for the halo distribution
double alfx
initial twiss parameters
double bety
initial twiss parameters
double haloPSWeightParameter
for the halo distribution
std::string haloPSWeightFunction
for the halo distribution
double haloXpCutInner
for the halo distribution
double haloYpCutOuter
for the halo distribution