BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSIntegratorDipoleQuadrupole.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 "BDSDebug.hh"
20#include "BDSFieldMagDipole.hh"
21#include "BDSIntegratorDipoleRodrigues2.hh"
22#include "BDSIntegratorDipoleQuadrupole.hh"
23#include "BDSIntegratorQuadrupole.hh"
24#include "BDSMagnetStrength.hh"
25#include "BDSParticleDefinition.hh"
26#include "BDSStep.hh"
27#include "BDSUtilities.hh"
28
29#include "globals.hh"
30#include "G4AffineTransform.hh"
31#include "G4Mag_EqRhs.hh"
32#include "G4MagIntegratorStepper.hh"
33#include "G4PhysicalConstants.hh"
34#include "G4SystemOfUnits.hh"
35#include "G4ThreeVector.hh"
36#include "G4Transform3D.hh"
37
38#include "CLHEP/Units/SystemOfUnits.h"
39
40#include <cmath>
41
43 G4double brhoIn,
44 G4Mag_EqRhs* eqOfMIn,
45 G4double minimumRadiusOfCurvatureIn,
46 const BDSParticleDefinition* designParticle,
47 const G4double& tiltIn):
48 BDSIntegratorMag(eqOfMIn, 6),
49 nominalBRho(brhoIn),
50 eq(dynamic_cast<BDSMagUsualEqRhs*>(eqOfM)),
51 bPrime(std::abs(brhoIn) * (*strengthIn)["k1"] / CLHEP::m2),
52 nominalBeta((*strengthIn)["beta0"]),
53 nominalRho((*strengthIn)["length"]/(*strengthIn)["angle"]),
54 nominalField((*strengthIn)["field"]),
55 fieldRatio(nominalField/ (nominalBRho/nominalRho)),
56 nominalEnergy(designParticle->TotalEnergy()),
57 nominalMass(designParticle->Mass()),
58 nominalCharge(designParticle->Charge()),
59 fieldArcLength((*strengthIn)["length"]),
60 nominalAngle((*strengthIn)["angle"]),
61 tilt(tiltIn),
62 scaling((*strengthIn)["scaling"]),
63 dipole(new BDSIntegratorDipoleRodrigues2(eqOfMIn, minimumRadiusOfCurvatureIn))
64{
65 if (!std::isfinite(nominalRho))
66 {nominalRho = 0;}
67 if (!std::isfinite(fieldRatio))
68 {fieldRatio = 1;}
69 isScaled = scaling == 1 ? false : true;
70 zeroStrength = !BDS::IsFinite((*strengthIn)["field"]);
71 BDSFieldMagDipole* dipoleField = new BDSFieldMagDipole(strengthIn);
72 unitField = (dipoleField->FieldValue()).unit();
73 delete dipoleField;
74 angleForCL = fieldRatio != 1 ? nominalAngle * fieldRatio : nominalAngle;
75 angleForCL /= scaling;
76 if (!std::isfinite(angleForCL))
77 {angleForCL = 0;}
78}
79
80BDSIntegratorDipoleQuadrupole::~BDSIntegratorDipoleQuadrupole()
81{
82 delete dipole;
83}
84
85void BDSIntegratorDipoleQuadrupole::Stepper(const G4double yIn[6],
86 const G4double dydx[6],
87 const G4double h,
88 G4double yOut[6],
89 G4double yErr[6])
90{
91 // charge and unit normalisation
92 const G4double fcof = eqOfM->FCof();
93
94 // Protect against very small steps, neutral particles, and zero field: drift through.
95 if (h < 1e-12 || !BDS::IsFiniteStrength(fcof) || zeroStrength)
96 {
97 AdvanceDriftMag(yIn,h,yOut,yErr);
98 SetDistChord(0);
99 return;
100 }
101
102 // try out a simple dipole step first - dumb stepper has no error but is 2x as fast as normal step
103 dipole->SingleStep(yIn, h, yOut);
104
105 G4double dipoleDC = dipole->DistChord();
106 // We assume that the major effect is the dipole component and the quadrupole
107 // component is low. Therefore, we can safely set distchord from the dipole.
108 SetDistChord(dipoleDC);
109
110 // return just dipole kick for small step size
111 if (h < 1e-4) // 1e-4mm
112 {
113 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
114 SetDistChord(dipole->DistChord());
115 return;
116 }
117
118 // If the particle might spiral, we return and just use the dipole only component
119 // Aimed at particles of much lower momentum than the design energy.
120 G4double radiusOfCurvature = dipole->RadiusOfHelix();
121 if (dipoleDC > 0.3*radiusOfCurvature)
122 {
123 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
124 SetDistChord(dipole->DistChord());
125 return;
126 }
127
128 // Revert to backup for zero angle, finite field dipoles. This could eventually be implemented with the
129 // matrix method by transforming onto the correct trajectory, but for high field, zero angle bends,
130 // the transform angle may be large. In this case, the distchord is large and the tracking is subsequently incorrect.
131 // As zero angle bends are constructed as one magnet segment, this should be solved by splitting the zero angle
132 // beampipe into segments like an ordinary sbend, but for now it's easier and safer to revert to backup.
134 {
135 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
136 SetDistChord(dipole->DistChord());
137 return;
138 }
139
140 // not going to spiral so proceed
141 // convert to true curvilinear
142 G4ThreeVector globalPos = G4ThreeVector(yIn[0], yIn[1], yIn[2]);
143 G4ThreeVector globalMom = G4ThreeVector(yIn[3], yIn[4], yIn[5]);
144
145 // create new position and momentum vectors for matrix output
146 G4ThreeVector localCLPosOut;
147 G4ThreeVector localCLMomOut;
148
149 // calculate relative energy difference
150 G4double nomMomentum = std::abs(nominalBRho * fcof); // safe as brho is nominal and abs(fcof) is scaling factor
151 G4double deltaEnergy = eq->TotalEnergy(globalMom) - nominalEnergy;
152 // deltaE/P0 to match literature.
153 G4double deltaEoverP0 = deltaEnergy / (nomMomentum);
154
155 // get particle charge from reverse of how it's calculated in G4Mag_EqRhs::SetChargeMomentumMass.
156 G4double pcharge = fcof/(eplus*c_light);
157 G4double chargeRatio = std::abs(pcharge /nominalCharge);
158
159 // revert to backup if off-charge
160 if (chargeRatio != 1.0)
161 {
162 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
163 SetDistChord(dipole->DistChord());
164 return;
165 }
166
167 if (isScaled)
168 {
169 // if the field is scaled to account for a central momentum change, the new nominal momentum must used
170 // in calculating the particle momentum offset. If the bunch's central momentum is not altered but the
171 // scaling is applied e.g for deflecting, then the particle will have a large relative momentum offset,
172 // therefore do not use the nominal matrix.
173 G4double nominalMomentum = std::abs(nominalBRho * fcof * fieldRatio);
174 G4double nomEnergy = std::sqrt(std::pow(nominalMomentum, 2) + std::pow(nominalMass,2));
175
176 // relative energy difference offset calculated w.r.t new nominal energy c/o scaling parameter
177 deltaEnergy = eq->TotalEnergy(globalMom) - nomEnergy;
178 // deltaE/P0 to match literature.
179 deltaEoverP0 = deltaEnergy / (nominalMomentum);
180 }
181
182 G4double beta = nominalBeta;
183 G4double rho = nominalRho;
184
185 if (isScaled)
186 {
187 G4double scaledEnergy = std::sqrt(std::pow(nomMomentum * scaling, 2) + std::pow(nominalMass, 2));
188 G4double scaledGamma = scaledEnergy / nominalMass;
189 beta = std::sqrt(1.0 - (1.0 / std::pow(scaledGamma, 2)));
190 }
191
192 // revert to backup if a particle is sufficiently off energy
193 if (std::abs(deltaEoverP0) > nominalMatrixRelativeMomCut)
194 {
195 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
196 SetDistChord(dipole->DistChord());
197 return;
198 //Leave this is as this may be re-investigated in the future.
199
200 /*
201 //transform onto non-nominal trajectory if particle is sufficiently off energy
202 //calculate angle for curvilinear transform
203 G4double particleBRho = globalMom.mag() / eq->FCof();
204 G4double scaleFactor = particleBRho / nominalBRho;
205 angleForTransform *= scaleFactor;
206
207 //calculate rho and beta for use in non-nominal matrix.
208 //rho according to field strength (from input) and particle momentum
209 rho = (globalMom.mag() / fcof) / nominalField;
210 beta = eq->Beta(yIn);
211 */
212 }
213
215 globalPos, globalMom, h, true, fcof, tilt);
216
217 G4ThreeVector localCLPos = localCL.PreStepPoint();
218 G4ThreeVector localCLMom = localCL.PostStepPoint();
219 G4ThreeVector localCLMomU = localCLMom.unit();
220
221 // only proceed with thick matrix if particle is paraxial
222 // judged by forward momentum > 1-limit and |transverse| < limit (default limit=0.1)
223 if (localCLMomU.z() < (1.0 - backupStepperMomLimit) || std::abs(localCLMomU.x()) > backupStepperMomLimit ||
224 std::abs(localCLMomU.y()) > backupStepperMomLimit)
225 {
226 dipole->Stepper(yIn, dydx, h, yOut, yErr); // more accurate step with error
227 SetDistChord(dipole->DistChord());
228 return;
229 }
230
231 // dipole step. Pass in rel. energy diff. to save recalculating.
232 OneStep(localCLPos, localCLMom, localCLMomU, h, fcof, localCLPosOut, localCLMomOut, rho, beta, deltaEoverP0);
233
234 // convert to global coordinates for output
235 BDSStep globalOut = CurvilinearToGlobal(fieldArcLength, unitField, angleForCL,
236 localCLPosOut, localCLMomOut, true, fcof, tilt);
237 G4ThreeVector globalPosOut = globalOut.PreStepPoint();
238 G4ThreeVector globalMomOut = globalOut.PostStepPoint();
239
240 // error along direction of travel really
241 G4ThreeVector globalMomOutU = globalMomOut.unit();
242 globalMomOutU *= 1e-10;
243
244 // chord distance (simple quadratic approx)
245 G4double dc = h*h/(8*radiusOfCurvature);
246 if (std::isnan(dc))
247 {SetDistChord(0);}
248 else
249 {SetDistChord(dc);}
250
251 // write out values and errors
252 for (G4int i = 0; i < 3; i++)
253 {
254 yOut[i] = globalPosOut[i];
255 yOut[i + 3] = globalMomOut[i];
256 yErr[i] = globalMomOutU[i];
257 yErr[i + 3] = 1e-40;
258 }
259}
260
261void BDSIntegratorDipoleQuadrupole::OneStep(const G4ThreeVector& posIn,
262 const G4ThreeVector& momIn,
263 const G4ThreeVector& momUIn,
264 const G4double& h,
265 const G4double& fcof,
266 G4ThreeVector& posOut,
267 G4ThreeVector& momOut,
268 const G4double rho,
269 const G4double beta,
270 const G4double deltaEoverP) const
271{
272 G4double momInMag = momIn.mag();
273 G4double nomMomentum = std::abs(nominalBRho * fcof); // safe as brho is nominal and abs(fcof) is scaling factor
274
275 // quad strength k normalised to charge and nominal momentum
276 // eqOfM->FCof() gives us conversion to MeV,mm and rigidity in Tm correctly
277 // as well as charge of the given particle
278 G4double K1 = std::abs(fcof)*bPrime/nomMomentum;
279 G4bool focussing = K1 >= 0; // depends on charge as well (in eqOfM->FCof())
280
281 // separate focussing strengths for vertical and horizontal axes.
282 // Used by matrix elements so must be derived from nominal values.
283 // strength parameter calculated differently depending on value
284 // default values are for +ve K1
285 G4double invrho2 = std::pow(1.0 / rho, 2);
286 G4double kx2 = invrho2 + K1;
287 G4double kx = std::sqrt(std::abs(kx2));
288 G4double kxl = kx * h;
289
290 G4double ky2 = -K1;
291 G4double ky = std::sqrt(std::abs(ky2));
292 G4double kyl = ky * h;
293
294 G4double x0 = posIn.x();
295 G4double y0 = posIn.y();
296 G4double s0 = posIn.z();
297 G4double xp = momUIn.x();
298 G4double yp = momUIn.y();
299 G4double zp = momUIn.z();
300
301 G4double x1 = x0;
302 G4double y1 = y0;
303 G4double xp1 = xp;
304 G4double yp1 = yp;
305 G4double zp1 = zp;
306
307 G4double X11=0,X12=0,X16=0,X21=0,X22=0,X26=0;
308 G4double Y11=0,Y12=0,Y21=0,Y22=0;
309
310 // matrix elements. All must derived from nominal parameters.
311 if (focussing)
312 {//focussing
313 X11 = std::cos(kxl);
314 X12 = std::sin(kxl)/kx;
315 X21 =-std::abs(kx2)*X12;
316 X22 = X11;
317 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cos(kxl));
318 X26 = (1.0/beta) * (1.0/rho) * X12;
319
320 Y11= std::cosh(kyl);
321 Y12= std::sinh(kyl)/ky;
322 if (std::isnan(Y12)) //y12 can be nan from div by 0 in previous line
323 {Y12 = h;}
324 Y21= std::abs(ky2)*Y12;
325 Y22= Y11;
326 }
327 else
328 {// defocussing
329 if (std::abs(K1) < invrho2)
330 {
331 kx2 = invrho2 - std::abs(K1);
332 kx = std::sqrt(std::abs(kx2));
333 kxl = kx * h;
334 X11 = std::cos(kxl);
335 X12 = std::sin(kxl)/kx; // always +ve
336 X21 = -std::abs(kx2)*X12;
337 X22 = X11;
338 }
339 else
340 {
341 kx2 = std::abs(K1) - invrho2;
342 kx = std::sqrt(std::abs(kx2));
343 kxl = kx * h;
344 X11 = std::cosh(kxl);
345 X12 = std::sinh(kxl)/kx; // always +ve
346 X21 = std::abs(kx2)*X12;
347 X22 = X11;
348 }
349 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cos(kxl));
350 X26 = (1.0/beta) * (1.0/rho) * X12;
351
352 Y11= std::cos(kyl);
353 Y12= std::sin(kyl)/ky;
354 if (std::isnan(Y12)) //y12 can be nan from div by 0 in previous line
355 {Y12 = h;}
356 Y21= -std::abs(ky2)*Y12;
357 Y22= Y11;
358 }
359
360 x1 = X11*x0 + X12*xp + X16*deltaEoverP;
361 xp1 = X21*x0 + X22*xp + X26*deltaEoverP;
362
363 y1 = Y11*y0 + Y12*yp;
364 yp1 = Y21*y0 + Y22*yp;
365
366 G4double s1 = s0 + h;
367
368 // relies on normalised momenta otherwise this will be nan.
369 zp1 = std::sqrt(1 - xp1*xp1 - yp1*yp1);
370 if (std::isnan(zp1))
371 {zp1 = zp;} // ensure not nan
372
373 G4ThreeVector momOutUnit = G4ThreeVector(xp1, yp1, zp1);
374 momOut = momOutUnit * momInMag;
375
376 posOut = G4ThreeVector(x1, y1, s1);
377
378}
379
BDSStep GlobalToCurvilinear(const G4double fieldArcLength, const G4ThreeVector &unitField, const G4double angle, const G4ThreeVector &position, const G4ThreeVector &unitMomentum, const G4double h, const G4bool useCurvilinearWorld, const G4double FCof, const G4double tilt=0)
A uniform dipole field.
G4double nominalRho
Cached magnet property, nominal bending radius.
G4double fieldRatio
Ratio of supplied field to nominal field. Needed for over/underpowered magnets.
virtual void Stepper(const G4double y[6], const G4double dydx[6], const G4double h, G4double yOut[6], G4double yErr[6])
const G4double fieldArcLength
Cache of the field arc length.
BDSIntegratorDipoleQuadrupole()=delete
Private default constructor to enforce use of supplied constructor.
void OneStep(const G4ThreeVector &posIn, const G4ThreeVector &momIn, const G4ThreeVector &momUIn, const G4double &h, const G4double &fcof, G4ThreeVector &posOut, G4ThreeVector &momOut, const G4double rho, const G4double beta, const G4double deltaEnergy) const
const G4double nominalCharge
Nominal beam charge.
const G4double bPrime
Cached magnet property, B field gradient for calculating K1.
const G4double nominalAngle
Cache of the field angle.
G4double angleForCL
Angle used for curvilinear transforms.
G4ThreeVector unitField
Cache of the unit field direction.
G4double tilt
Tilt offset transform for field.
const G4double nominalMass
Primary particle mass. Needed for recalculating nominal energy with scaling.
const G4double nominalBeta
Cached nominal relativistic beta of the nominal beam particle.
const G4double scaling
Cache field scaling factor.
G4bool isScaled
Cache if field is scaled.
BDSMagUsualEqRhs * eq
BDSIM's eqRhs class to give access to particle properties.
const G4double nominalBRho
Cached magnet property, nominal magnetic rigidity.
const G4double nominalEnergy
Nominal beam energy.
BDSIntegratorDipoleRodrigues2 * dipole
Backup integrator.
Exact helix through pure dipole field.
void SingleStep(const G4double yIn[6], const G4double &h, G4double yOut[6])
G4double RadiusOfHelix() const
Public accessor for protected variable in base class.
virtual void Stepper(const G4double yIn[6], const G4double dydx[], G4double h, G4double yOut[6], G4double yErr[6])
Output error array.
void AdvanceDriftMag(const G4double yIn[], const G4double h, G4double yOut[], G4double yErr[]) const
Error array [6] all 0.
Common functionality to BDSIM integrators.
G4double backupStepperMomLimit
static G4double nominalMatrixRelativeMomCut
G4Mag_EqRhs * eqOfM
void SetDistChord(G4double distChordIn)
Setter for distChord to private member.
Override G4Mag_UsualEqRhs, provides BDSIM integrators access to particle attributes.
G4double TotalEnergy(const G4double y[])
Calculate total particle energy.
Efficient storage of magnet strengths.
Wrapper for particle definition.
A simple class to represent the positions of a step.
Definition: BDSStep.hh:33
G4ThreeVector PostStepPoint() const
Accessor.
Definition: BDSStep.hh:43
G4ThreeVector PreStepPoint() const
Accessor.
Definition: BDSStep.hh:42
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())
G4bool IsFiniteStrength(G4double variable)