20#include "BDSFieldMagDipole.hh"
21#include "BDSIntegratorDipoleRodrigues2.hh"
22#include "BDSIntegratorDipoleQuadrupole.hh"
23#include "BDSIntegratorQuadrupole.hh"
24#include "BDSMagnetStrength.hh"
25#include "BDSParticleDefinition.hh"
27#include "BDSUtilities.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"
38#include "CLHEP/Units/SystemOfUnits.h"
45 G4double minimumRadiusOfCurvatureIn,
47 const G4double& tiltIn):
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"]),
62 scaling((*strengthIn)[
"scaling"]),
65 if (!std::isfinite(nominalRho))
67 if (!std::isfinite(fieldRatio))
69 isScaled = scaling == 1 ? false :
true;
72 unitField = (dipoleField->FieldValue()).unit();
74 angleForCL = fieldRatio != 1 ? nominalAngle * fieldRatio : nominalAngle;
75 angleForCL /= scaling;
76 if (!std::isfinite(angleForCL))
80BDSIntegratorDipoleQuadrupole::~BDSIntegratorDipoleQuadrupole()
86 const G4double dydx[6],
92 const G4double fcof =
eqOfM->FCof();
105 G4double dipoleDC =
dipole->DistChord();
121 if (dipoleDC > 0.3*radiusOfCurvature)
142 G4ThreeVector globalPos = G4ThreeVector(yIn[0], yIn[1], yIn[2]);
143 G4ThreeVector globalMom = G4ThreeVector(yIn[3], yIn[4], yIn[5]);
146 G4ThreeVector localCLPosOut;
147 G4ThreeVector localCLMomOut;
150 G4double nomMomentum = std::abs(
nominalBRho * fcof);
153 G4double deltaEoverP0 = deltaEnergy / (nomMomentum);
156 G4double pcharge = fcof/(eplus*c_light);
160 if (chargeRatio != 1.0)
174 G4double nomEnergy = std::sqrt(std::pow(nominalMomentum, 2) + std::pow(
nominalMass,2));
179 deltaEoverP0 = deltaEnergy / (nominalMomentum);
187 G4double scaledEnergy = std::sqrt(std::pow(nomMomentum *
scaling, 2) + std::pow(
nominalMass, 2));
189 beta = std::sqrt(1.0 - (1.0 / std::pow(scaledGamma, 2)));
215 globalPos, globalMom, h,
true, fcof,
tilt);
219 G4ThreeVector localCLMomU = localCLMom.unit();
232 OneStep(localCLPos, localCLMom, localCLMomU, h, fcof, localCLPosOut, localCLMomOut, rho, beta, deltaEoverP0);
236 localCLPosOut, localCLMomOut,
true, fcof,
tilt);
241 G4ThreeVector globalMomOutU = globalMomOut.unit();
242 globalMomOutU *= 1e-10;
245 G4double dc = h*h/(8*radiusOfCurvature);
252 for (G4int i = 0; i < 3; i++)
254 yOut[i] = globalPosOut[i];
255 yOut[i + 3] = globalMomOut[i];
256 yErr[i] = globalMomOutU[i];
262 const G4ThreeVector& momIn,
263 const G4ThreeVector& momUIn,
265 const G4double& fcof,
266 G4ThreeVector& posOut,
267 G4ThreeVector& momOut,
270 const G4double deltaEoverP)
const
272 G4double momInMag = momIn.mag();
273 G4double nomMomentum = std::abs(
nominalBRho * fcof);
278 G4double K1 = std::abs(fcof)*
bPrime/nomMomentum;
279 G4bool focussing = K1 >= 0;
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;
291 G4double ky = std::sqrt(std::abs(ky2));
292 G4double kyl = ky * h;
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();
307 G4double X11=0,X12=0,X16=0,X21=0,X22=0,X26=0;
308 G4double Y11=0,Y12=0,Y21=0,Y22=0;
314 X12 = std::sin(kxl)/kx;
315 X21 =-std::abs(kx2)*X12;
317 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cos(kxl));
318 X26 = (1.0/beta) * (1.0/rho) * X12;
321 Y12= std::sinh(kyl)/ky;
324 Y21= std::abs(ky2)*Y12;
329 if (std::abs(K1) < invrho2)
331 kx2 = invrho2 - std::abs(K1);
332 kx = std::sqrt(std::abs(kx2));
335 X12 = std::sin(kxl)/kx;
336 X21 = -std::abs(kx2)*X12;
341 kx2 = std::abs(K1) - invrho2;
342 kx = std::sqrt(std::abs(kx2));
344 X11 = std::cosh(kxl);
345 X12 = std::sinh(kxl)/kx;
346 X21 = std::abs(kx2)*X12;
349 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cos(kxl));
350 X26 = (1.0/beta) * (1.0/rho) * X12;
353 Y12= std::sin(kyl)/ky;
356 Y21= -std::abs(ky2)*Y12;
360 x1 = X11*x0 + X12*xp + X16*deltaEoverP;
361 xp1 = X21*x0 + X22*xp + X26*deltaEoverP;
363 y1 = Y11*y0 + Y12*yp;
364 yp1 = Y21*y0 + Y22*yp;
366 G4double s1 = s0 + h;
369 zp1 = std::sqrt(1 - xp1*xp1 - yp1*yp1);
373 G4ThreeVector momOutUnit = G4ThreeVector(xp1, yp1, zp1);
374 momOut = momOutUnit * momInMag;
376 posOut = G4ThreeVector(x1, y1, s1);
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)
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
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.
G4ThreeVector PostStepPoint() const
Accessor.
G4ThreeVector PreStepPoint() const
Accessor.
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())
G4bool IsFiniteStrength(G4double variable)