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)
60BDSBunchHalo::~BDSBunchHalo()
66 G4Transform3D beamlineTransformIn,
67 const G4double beamlineSIn)
99 emitInnerX = std::pow(haloNSigmaXInner, 2) *
emitX;
100 emitInnerY = std::pow(haloNSigmaYInner, 2) *
emitY;
101 emitOuterX = std::pow(haloNSigmaXOuter, 2) *
emitX;
102 emitOuterY = std::pow(haloNSigmaYOuter, 2) *
emitY;
104 xMax = haloNSigmaXOuter *
sigmaX;
105 yMax = haloNSigmaYOuter *
sigmaY;
106 xpMax = std::sqrt(std::pow(haloNSigmaXOuter, 2) *
emitX *
gammaX);
107 ypMax = std::sqrt(std::pow(haloNSigmaYOuter, 2) *
emitY *
gammaY);
126 G4double dx = xMax * (1 - 2 * G4RandFlat::shoot());
127 G4double dy = yMax * (1 - 2 * G4RandFlat::shoot());
128 G4double dxp = xpMax * (1 - 2 * G4RandFlat::shoot());
129 G4double dyp = ypMax * (1 - 2 * G4RandFlat::shoot());
132 double emitXSp =
gammaX * std::pow(std::abs(dx), 2) + (2. *
alphaX * dx * dxp) +
betaX * std::pow(std::abs(dxp), 2);
133 double emitYSp =
gammaY * std::pow(std::abs(dy), 2) + (2. *
alphaY * dy * dyp) +
betaY * std::pow(std::abs(dyp), 2);
137 if ((std::abs(emitXSp) < emitInnerX || std::abs(emitYSp) < emitInnerY) ||
138 (std::abs(emitXSp) > emitOuterX || std::abs(emitYSp) > emitOuterY) ||
139 (std::abs(dx) < (haloXCutInner *
sigmaX)) ||
140 (std::abs(dy) < (haloYCutInner *
sigmaY)) ||
141 (std::abs(dx) > (haloXCutOuter *
sigmaX)) ||
142 (std::abs(dy) > (haloYCutOuter *
sigmaY)) ||
143 (std::abs(dxp) < (haloXpCutInner *
sigmaXp)) ||
144 (std::abs(dyp) < (haloYpCutInner *
sigmaYp)) ||
145 (std::abs(dxp) > (haloXpCutOuter *
sigmaXp)) ||
146 (std::abs(dyp) > (haloYpCutOuter *
sigmaYp)) )
155 if (weightFunction ==
"flat" || weightFunction.empty() || weightFunction ==
"one")
160 else if (weightFunction ==
"oneoverr")
163 wx = std::pow(std::abs(emitInnerX / emitXSp), haloPSWeightParameter);
164 wy = std::pow(std::abs(emitInnerY / emitYSp), haloPSWeightParameter);
166 else if (weightFunction ==
"oneoverrsqrd")
169 double eXsqrd = std::pow(std::abs(emitXSp), 2);
170 double eYsqrd = std::pow(std::abs(emitYSp), 2);
171 double eXInsq = std::pow(std::abs(emitInnerX), 2);
172 double eYInsq = std::pow(std::abs(emitInnerY), 2);
173 wx = std::pow(std::abs(eXInsq / eXsqrd), haloPSWeightParameter);
174 wy = std::pow(std::abs(eYInsq / eYsqrd), haloPSWeightParameter);
176 else if (weightFunction ==
"exp")
178 wx = std::exp(-(emitXSp * haloPSWeightParameter) / (emitInnerX));
179 wy = std::exp(-(emitYSp * haloPSWeightParameter) / (emitInnerY));
183 G4cout << __METHOD_NAME__ << emitXSp/
emitX <<
" " << emitYSp/
emitY <<
" " << wx <<
" " << wy << G4endl;
186 if (G4RandFlat::shoot() > wx && G4RandFlat::shoot() > wy)
192 xp += dxp * CLHEP::rad;
193 yp += dyp * CLHEP::rad;
198 G4cout << __METHOD_NAME__ <<
"selected> " << dx <<
" " << dy <<
" " << dxp <<
" " << dyp << G4endl;
212 {
throw BDSException(__METHOD_NAME__,
"emitx must be finite!");}
214 {
throw BDSException(__METHOD_NAME__,
"emity must be finite!");}
217 {
throw BDSException(__METHOD_NAME__,
"betx must be finite!");}
219 {
throw BDSException(__METHOD_NAME__,
"bety must be finite!");}
221 std::vector<G4String> weightFunctions = {
"",
"one",
"flat",
"oneoverr",
"oneoverrsqrd",
"exp"};
222 auto search = std::find(weightFunctions.begin(), weightFunctions.end(), weightFunction);
223 if (search == weightFunctions.end())
225 G4cerr << __METHOD_END__ <<
"invalid haloPSWeightFunction \"" << weightFunction <<
"\"" << G4endl;
226 G4cout <<
"Available weight functions are:" << G4endl;
227 for (
const auto& w : weightFunctions)
228 {G4cout << w << G4endl;}
232 if (haloNSigmaXInner <= 0)
233 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner <= 0");}
235 if (haloNSigmaYInner <= 0)
236 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner <= 0");}
238 if (haloNSigmaXInner > haloNSigmaXOuter)
239 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaXInner cannot be less than haloNSigmaXOuter");}
241 if (haloNSigmaYInner > haloNSigmaYOuter)
242 {
throw BDSException(__METHOD_NAME__,
"haloNSigmaYInner cannot be less than haloNSigmaYOuter");}
244 if (haloXCutInner < 0)
245 {
throw BDSException(__METHOD_NAME__,
"haloXCutInner < 0");}
247 if (haloYCutInner < 0)
248 {
throw BDSException(__METHOD_NAME__,
"haloYCutInner < 0");}
250 if (haloXCutOuter <= haloXCutInner)
251 {
throw BDSException(__METHOD_NAME__,
"haloXCutOuter must be greater than haloXCutInner!");}
253 if (haloYCutOuter <= haloYCutInner)
254 {
throw BDSException(__METHOD_NAME__,
"haloYCutOuter must be greater than haloYCutInner!");}
256 if (haloXpCutInner < 0)
257 {
throw BDSException(__METHOD_NAME__,
"haloXpCutInner < 0");}
259 if (haloYpCutInner < 0)
260 {
throw BDSException(__METHOD_NAME__,
"haloYpCutInner < 0");}
262 if (haloXpCutOuter <= haloXpCutInner)
263 {
throw BDSException(__METHOD_NAME__,
"haloXpCutOuter must be greater than haloXpCutInner!");}
265 if (haloYpCutOuter <= haloYpCutInner)
266 {
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