BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSIntegratorSextupole.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 "BDSIntegratorSextupole.hh"
21#include "BDSMagnetStrength.hh"
22#include "BDSStep.hh"
23#include "BDSUtilities.hh"
24
25#include "globals.hh" // geant4 types / globals
26#include "G4Mag_EqRhs.hh"
27#include "G4ThreeVector.hh"
28
29#include "CLHEP/Units/SystemOfUnits.h"
30
31#include <cmath>
32
33BDSIntegratorSextupole::BDSIntegratorSextupole(BDSMagnetStrength const* strength,
34 G4double brho,
35 G4Mag_EqRhs* eqOfMIn):
37{
38 // B'' = d^2By/dx^2 = Brho * (1/Brho d^2By/dx^2) = Brho * k2
39 bDoublePrime = brho * (*strength)["k2"] / CLHEP::m3;
40 zeroStrength = !BDS::IsFinite(bDoublePrime);
41#ifdef BDSDEBUG
42 G4cout << __METHOD_NAME__ << "B'' = " << bDoublePrime << G4endl;
43#endif
44}
45
46void BDSIntegratorSextupole::AdvanceHelix(const G4double yIn[],
47 G4double h,
48 G4double yOut[],
49 G4double yErr[])
50{
51 G4ThreeVector mom = G4ThreeVector(yIn[3], yIn[4], yIn[5]);
52 G4double momUnit = mom.mag();
53 G4double kappa = (-eqOfM->FCof()*bDoublePrime) / momUnit;
54
55 if (std::abs(kappa) < 1e-12)
56 {
57 AdvanceDriftMag(yIn, h, yOut, yErr);
58 SetDistChord(0);
59 return;
60 }
61
62 G4ThreeVector pos = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
63 BDSStep localPosMom = ConvertToLocal(pos, mom, h, false);
64 G4ThreeVector localPos = localPosMom.PreStepPoint();
65 G4ThreeVector localMom = localPosMom.PostStepPoint();
66 G4ThreeVector localMomUnit = localMom.unit();
67
68 G4double x0 = localPos.x();
69 G4double y0 = localPos.y();
70
71 G4double xp = localMomUnit.x();
72 G4double yp = localMomUnit.y();
73 G4double zp = localMomUnit.z();
74
75 // Evaluate field at the approximate midpoint of the step.
76 const G4double halfH = 0.5*h;
77 x0 = x0 + localMomUnit.x()*halfH;
78 y0 = y0 + localMomUnit.y()*halfH;
79
80 G4double x02My02 = x0*x0 - y0*y0;
81
82 // local r'' (for curvature)
83 G4ThreeVector localA;
84 localA.setX(zp*x02My02);
85 localA.setY(-2*zp*x0*y0);
86 localA.setZ(xp*x02My02-2*yp*x0*y0);
87
88 localA *= kappa / 2; // 2 is actually a 2! factor.
89
90 AdvanceChord(h,localPos,localMomUnit,localA);
91 ConvertToGlobal(localPos, localMomUnit, yOut, yErr, momUnit);
92}
BDSStep ConvertToLocal(G4Step const *const step, G4bool useCurvilinear=true) const
void AdvanceDriftMag(const G4double yIn[], const G4double h, G4double yOut[], G4double yErr[]) const
Error array [6] all 0.
Common functionality for Euler integrators.
void AdvanceChord(const G4double h, G4ThreeVector &localPos, G4ThreeVector &localMom, const G4ThreeVector &localA)
G4Mag_EqRhs * eqOfM
void SetDistChord(G4double distChordIn)
Setter for distChord to private member.
void ConvertToGlobal(const G4ThreeVector &localPos, const G4ThreeVector &localMom, G4double yOut[], G4double yErr[], const G4double momScaling=1.0)
scaling of momentum in case localMom is a unit vector
virtual void AdvanceHelix(const G4double yIn[], G4double h, G4double yOut[], G4double yErr[])
Calculate the new particle coordinates.
G4double bDoublePrime
2nd derivative of the field
Efficient storage of magnet strengths.
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())