BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSBunchHalo.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSBunchHalo.hh"
20#include "BDSDebug.hh"
21#include "BDSException.hh"
22#include "BDSParticleCoordsFull.hh"
23
24#include "parser/beam.h"
25
26#include "globals.hh" // geant4 types / globals
27
28#include "Randomize.hh"
29#include "CLHEP/Random/RandFlat.h"
30#include "CLHEP/Units/SystemOfUnits.h"
31
32#include <cmath>
33
34BDSBunchHalo::BDSBunchHalo():
35 BDSBunch("halo"),
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),
44 haloXCutInner(0.0),
45 haloYCutInner(0.0),
46 haloXCutOuter(0.0),
47 haloYCutOuter(0.0),
48 haloXpCutInner(0.0),
49 haloYpCutInner(0.0),
50 haloXpCutOuter(0.0),
51 haloYpCutOuter(0.0),
52 haloPSWeightParameter(0.0),
53 weightFunction(""),
54 emitInnerX(0.0), emitInnerY(0.0),
55 emitOuterX(0.0), emitOuterY(0.0),
56 xMax(0.0), yMax(0.0),
57 xpMax(0.0), ypMax(0.0)
58{;}
59
60BDSBunchHalo::~BDSBunchHalo()
61{;}
62
64 const GMAD::Beam& beam,
65 const BDSBunchType& distrType,
66 G4Transform3D beamlineTransformIn,
67 const G4double beamlineSIn)
68{
69 BDSBunch::SetOptions(beamParticle, beam, distrType, beamlineTransformIn, beamlineSIn);
70 alphaX = G4double(beam.alfx);
71 alphaY = G4double(beam.alfy);
72 betaX = G4double(beam.betx);
73 betaY = G4double(beam.bety);
74 gammaX = (1.0+alphaX*alphaX)/betaX;
75 gammaY = (1.0+alphaY*alphaY)/betaY;
76 haloNSigmaXInner = G4double(beam.haloNSigmaXInner);
77 haloNSigmaXOuter = G4double(beam.haloNSigmaXOuter);
78 haloNSigmaYInner = G4double(beam.haloNSigmaYInner);
79 haloNSigmaYOuter = G4double(beam.haloNSigmaYOuter);
80 haloXCutInner = G4double(beam.haloXCutInner);
81 haloYCutInner = G4double(beam.haloYCutInner);
82 haloXCutOuter = G4double(beam.haloXCutOuter);
83 haloYCutOuter = G4double(beam.haloYCutOuter);
84 haloXpCutInner = G4double(beam.haloXpCutInner);
85 haloYpCutInner = G4double(beam.haloYpCutInner);
86 haloXpCutOuter = G4double(beam.haloXpCutOuter);
87 haloYpCutOuter = G4double(beam.haloYpCutOuter);
88 haloPSWeightParameter = G4double(beam.haloPSWeightParameter);
89 weightFunction = G4String(beam.haloPSWeightFunction);
90
91 G4double ex,ey; // dummy variables we don't need
92 SetEmittances(beamParticle, beam, emitX, emitY, ex, ey);
93
94 sigmaX = std::sqrt(emitX * betaX);
95 sigmaY = std::sqrt(emitY * betaY);
96 sigmaXp = std::sqrt(gammaX * emitX);
97 sigmaYp = std::sqrt(gammaY * emitY);
98
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;
103
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);
108
110}
111
113{
114 // Central orbit
115 G4double x = X0;
116 G4double y = Y0;
117 G4double z = Z0;
118 G4double xp = Xp0;
119 G4double yp = Yp0;
120
121 // z += (T0 - envelopeT * (1.-2.*flatGen->shoot())) * CLHEP::c_light * CLHEP::s;
122 z = 0;
123
124 while (true)
125 {
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());
130
131 // compute single particle emittance
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);
134
135 // check if particle is within normal beam core, if so continue generation
136 // also check if particle is within the desired cut.
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)) )
147 {
148 continue;
149 }
150 else
151 {
152 // determine weight, initialise 1 so always passes
153 double wx = 1.0;
154 double wy = 1.0;
155 if (weightFunction == "flat" || weightFunction.empty() || weightFunction == "one")
156 {
157 wx = 1.0;
158 wy = 1.0;
159 }
160 else if (weightFunction == "oneoverr")
161 {
162 //abs because power of double - must be positive
163 wx = std::pow(std::abs(emitInnerX / emitXSp), haloPSWeightParameter);
164 wy = std::pow(std::abs(emitInnerY / emitYSp), haloPSWeightParameter);
165 }
166 else if (weightFunction == "oneoverrsqrd")
167 {
168 //abs because power of double - must be positive
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);
175 }
176 else if (weightFunction == "exp")
177 {
178 wx = std::exp(-(emitXSp * haloPSWeightParameter) / (emitInnerX));
179 wy = std::exp(-(emitYSp * haloPSWeightParameter) / (emitInnerY));
180 }
181
182#ifdef BDSDEBUG
183 G4cout << __METHOD_NAME__ << emitXSp/emitX << " " << emitYSp/emitY << " " << wx << " " << wy << G4endl;
184#endif
185 // reject
186 if (G4RandFlat::shoot() > wx && G4RandFlat::shoot() > wy)
187 {continue;}
188
189 // add to reference orbit
190 x += dx * CLHEP::m;
191 y += dy * CLHEP::m;
192 xp += dxp * CLHEP::rad;
193 yp += dyp * CLHEP::rad;
194
195 G4double zp = CalculateZp(xp, yp, Zp0);
196
197#ifdef BDSDEBUG
198 G4cout << __METHOD_NAME__ << "selected> " << dx << " " << dy << " " << dxp << " " << dyp << G4endl;
199#endif
200 // E0 and T0 from base class and already in G4 units
201 BDSParticleCoordsFull result(x,y,z,xp,yp,zp,T0,S0+z,E0,/*weight=*/1.0);
202 return result;
203 }
204 }
205}
206
208{
210
211 if (emitX <= 0)
212 {throw BDSException(__METHOD_NAME__, "emitx must be finite!");}
213 if (emitY <= 0)
214 {throw BDSException(__METHOD_NAME__, "emity must be finite!");}
215
216 if (betaX <= 0)
217 {throw BDSException(__METHOD_NAME__, "betx must be finite!");}
218 if (betaY <= 0)
219 {throw BDSException(__METHOD_NAME__, "bety must be finite!");}
220
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())
224 {
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;}
229 throw BDSException(__METHOD_END__, "");
230 }
231
232 if (haloNSigmaXInner <= 0)
233 {throw BDSException(__METHOD_NAME__, "haloNSigmaXInner <= 0");}
234
235 if (haloNSigmaYInner <= 0)
236 {throw BDSException(__METHOD_NAME__, "haloNSigmaYInner <= 0");}
237
238 if (haloNSigmaXInner > haloNSigmaXOuter)
239 {throw BDSException(__METHOD_NAME__, "haloNSigmaXInner cannot be less than haloNSigmaXOuter");}
240
241 if (haloNSigmaYInner > haloNSigmaYOuter)
242 {throw BDSException(__METHOD_NAME__, "haloNSigmaYInner cannot be less than haloNSigmaYOuter");}
243
244 if (haloXCutInner < 0)
245 {throw BDSException(__METHOD_NAME__, "haloXCutInner < 0");}
246
247 if (haloYCutInner < 0)
248 {throw BDSException(__METHOD_NAME__, "haloYCutInner < 0");}
249
250 if (haloXCutOuter <= haloXCutInner)
251 {throw BDSException(__METHOD_NAME__, "haloXCutOuter must be greater than haloXCutInner!");}
252
253 if (haloYCutOuter <= haloYCutInner)
254 {throw BDSException(__METHOD_NAME__, "haloYCutOuter must be greater than haloYCutInner!");}
255
256 if (haloXpCutInner < 0)
257 {throw BDSException(__METHOD_NAME__, "haloXpCutInner < 0");}
258
259 if (haloYpCutInner < 0)
260 {throw BDSException(__METHOD_NAME__, "haloYpCutInner < 0");}
261
262 if (haloXpCutOuter <= haloXpCutInner)
263 {throw BDSException(__METHOD_NAME__, "haloXpCutOuter must be greater than haloXpCutInner!");}
264
265 if (haloYpCutOuter <= haloYpCutInner)
266 {throw BDSException(__METHOD_NAME__, "haloYpCutOuter must be greater than haloYpCutInner!");}
267}
G4double sigmaX
Twiss parameter.
Definition: BDSBunchHalo.hh:61
G4double sigmaYp
Twiss parameter.
Definition: BDSBunchHalo.hh:64
G4double alphaY
Twiss parameter.
Definition: BDSBunchHalo.hh:54
G4double gammaY
Twiss parameter.
Definition: BDSBunchHalo.hh:60
G4double emitX
Twiss parameter.
Definition: BDSBunchHalo.hh:57
G4double gammaX
Twiss parameter.
Definition: BDSBunchHalo.hh:59
G4double sigmaXp
Twiss parameter.
Definition: BDSBunchHalo.hh:63
G4double emitY
Twiss parameter.
Definition: BDSBunchHalo.hh:58
G4double sigmaY
Twiss parameter.
Definition: BDSBunchHalo.hh:62
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
Definition: BDSBunchHalo.cc:63
virtual BDSParticleCoordsFull GetNextParticleLocal()
G4double betaY
Twiss parameter.
Definition: BDSBunchHalo.hh:56
virtual void CheckParameters()
G4double alphaX
Twiss parameter.
Definition: BDSBunchHalo.hh:53
G4double betaX
Twiss parameter.
Definition: BDSBunchHalo.hh:55
The base class for bunch distribution generators.
Definition: BDSBunch.hh:47
G4double Yp0
Centre of distributions.
Definition: BDSBunch.hh:177
G4double T0
Centre of distributions.
Definition: BDSBunch.hh:175
G4double S0
Centre of distributions.
Definition: BDSBunch.hh:174
G4double Z0
Centre of distributions.
Definition: BDSBunch.hh:173
G4double X0
Centre of distributions.
Definition: BDSBunch.hh:171
G4double Zp0
Centre of distributions.
Definition: BDSBunch.hh:178
G4double Xp0
Centre of distributions.
Definition: BDSBunch.hh:176
G4double E0
Centre of distributions.
Definition: BDSBunch.hh:179
static void SetEmittances(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, G4double &emittGeometricX, G4double &emittGeometricY, G4double &emittNormalisedX, G4double &emittNormalisedY)
Definition: BDSBunch.cc:175
G4double Y0
Centre of distributions.
Definition: BDSBunch.hh:172
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
Definition: BDSBunch.cc:82
static G4double CalculateZp(G4double xp, G4double yp, G4double Zp0)
Calculate zp safely based on other components.
Definition: BDSBunch.cc:367
virtual void CheckParameters()
Definition: BDSBunch.cc:223
General exception with possible name of object and message.
Definition: BDSException.hh:35
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
Definition: beamBase.h:90
double haloXCutOuter
for the halo distribution
Definition: beamBase.h:130
double haloNSigmaYInner
for the halo distribution
Definition: beamBase.h:126
double haloXpCutOuter
for the halo distribution
Definition: beamBase.h:134
double haloXCutInner
for the halo distribution
Definition: beamBase.h:128
double haloYpCutInner
for the halo distribution
Definition: beamBase.h:133
double betx
initial twiss parameters
Definition: beamBase.h:90
double haloYCutOuter
for the halo distribution
Definition: beamBase.h:131
double haloNSigmaXInner
for the halo distribution
Definition: beamBase.h:124
double haloNSigmaYOuter
for the halo distribution
Definition: beamBase.h:127
double haloYCutInner
for the halo distribution
Definition: beamBase.h:129
double haloNSigmaXOuter
for the halo distribution
Definition: beamBase.h:125
double alfx
initial twiss parameters
Definition: beamBase.h:90
double bety
initial twiss parameters
Definition: beamBase.h:90
double haloPSWeightParameter
for the halo distribution
Definition: beamBase.h:136
std::string haloPSWeightFunction
for the halo distribution
Definition: beamBase.h:137
double haloXpCutInner
for the halo distribution
Definition: beamBase.h:132
double haloYpCutOuter
for the halo distribution
Definition: beamBase.h:135
Beam class.
Definition: beam.h:44