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;
282 G4double kx2 = std::pow(1.0 / rho, 2) + K1;
283 G4double kx = std::sqrt(std::abs(kx2));
285 G4double ky = std::sqrt(std::abs(ky2));
286 G4double kxl = kx * h;
287 G4double kyl = ky * h;
289 G4bool focussing = K1 >= 0;
291 G4double x0 = posIn.x();
292 G4double y0 = posIn.y();
293 G4double s0 = posIn.z();
294 G4double xp = momUIn.x();
295 G4double yp = momUIn.y();
296 G4double zp = momUIn.z();
304 G4double X11=0,X12=0,X16=0,X21=0,X22=0,X26=0;
305 G4double Y11=0,Y12=0,Y21=0,Y22=0;
311 X12= std::sin(kxl)/kx;
312 X21=-std::abs(kx2)*X12;
314 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cos(kxl));
315 X26 = (1.0/beta) * (1.0/rho) * X12;
318 Y12= std::sinh(kyl)/ky;
321 Y21= std::abs(ky2)*Y12;
327 X12= std::sinh(kxl)/kx;
328 X21= std::abs(kx2)*X12;
330 X16 = (1.0/beta) * ((1.0/rho) / kx2) * (1 - std::cosh(kxl));
331 X26 = (1.0/beta) * (1.0/rho) * X12;
334 Y12= std::sin(kyl)/ky;
337 Y21= -std::abs(ky2)*Y12;
341 x1 = X11*x0 + X12*xp + X16*deltaEoverP;
342 xp1 = X21*x0 + X22*xp + X26*deltaEoverP;
344 y1 = Y11*y0 + Y12*yp;
345 yp1 = Y21*y0 + Y22*yp;
347 G4double s1 = s0 + h;
350 zp1 = std::sqrt(1 - xp1*xp1 - yp1*yp1);
354 G4ThreeVector momOutUnit = G4ThreeVector(xp1, yp1, zp1);
355 momOut = momOutUnit * momInMag;
357 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)