BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSFieldInfo.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 "BDSArrayReflectionType.hh"
20#include "BDSFieldInfo.hh"
21#include "BDSFieldType.hh"
22#include "BDSIntegratorType.hh"
23#include "BDSInterpolatorType.hh"
24#include "BDSMagnetStrength.hh"
25#include "BDSModulatorInfo.hh"
26#include "BDSUtilities.hh"
27
28#include "globals.hh" // geant4 types / globals
29#include "G4RotationMatrix.hh"
30#include "G4ThreeVector.hh"
31#include "G4Track.hh"
32#include "G4Transform3D.hh"
33#include "G4UserLimits.hh"
34
35#include <algorithm>
36#include <ostream>
37
38G4UserLimits* BDSFieldInfo::defaultUL = nullptr;
39
41 fieldType(BDSFieldType::none),
42 brho(0),
43 integratorType(BDSIntegratorType::none),
44 magnetStrength(nullptr),
45 provideGlobalTransform(false),
46 transform(nullptr),
47 magneticFieldFilePath(""),
48 magneticFieldFormat(BDSFieldFormat::none),
49 magneticInterpolatorType(BDSInterpolatorType::nearest3d),
50 magneticArrayReflectionTypeSet(BDSArrayReflectionTypeSet()),
51 electricFieldFilePath(""),
52 electricFieldFormat(BDSFieldFormat::none),
53 electricInterpolatorType(BDSInterpolatorType::nearest3d),
54 electricArrayReflectionTypeSet(BDSArrayReflectionTypeSet()),
55 cacheTransforms(true),
56 eScaling(1.0),
57 bScaling(1.0),
58 timeOffset(0),
59 autoScale(false),
60 stepLimit(nullptr),
61 poleTipRadius(1),
62 beamPipeRadius(0),
63 chordStepMinimum(-1),
64 tilt(0),
65 secondFieldOnLeft(false),
66 magneticSubFieldName(""),
67 electricSubFieldName(""),
68 usePlacementWorldTransform(false),
69 modulatorInfo(nullptr),
70 transformBeamline(nullptr),
71 nameOfParserDefinition("")
72{;}
73
75 G4double brhoIn,
76 BDSIntegratorType integratorTypeIn,
77 BDSMagnetStrength* magnetStrengthIn,
78 G4bool provideGlobalTransformIn,
79 const G4Transform3D& transformIn,
80 const G4String& magneticFieldFilePathIn,
81 BDSFieldFormat magneticFieldFormatIn,
82 BDSInterpolatorType magneticInterpolatorTypeIn,
83 const G4String& electricFieldFilePathIn,
84 BDSFieldFormat electricFieldFormatIn,
85 BDSInterpolatorType electricInterpolatorTypeIn,
86 G4bool cacheTransformsIn,
87 G4double eScalingIn,
88 G4double bScalingIn,
89 G4double timeOffsetIn,
90 G4bool autoScaleIn,
91 G4UserLimits* stepLimitIn,
92 G4double poleTipRadiusIn,
93 G4double beamPipeRadiusIn,
94 G4bool leftIn,
95 const G4String& magneticSubFieldNameIn,
96 const G4String& electricSubFieldNameIn):
97 fieldType(fieldTypeIn),
98 brho(brhoIn),
99 integratorType(integratorTypeIn),
100 magnetStrength(magnetStrengthIn),
101 provideGlobalTransform(provideGlobalTransformIn),
102 transform(nullptr),
103 magneticFieldFilePath(magneticFieldFilePathIn),
104 magneticFieldFormat(magneticFieldFormatIn),
105 magneticInterpolatorType(magneticInterpolatorTypeIn),
106 magneticArrayReflectionTypeSet(BDSArrayReflectionTypeSet()),
107 electricFieldFilePath(electricFieldFilePathIn),
108 electricFieldFormat(electricFieldFormatIn),
109 electricInterpolatorType(electricInterpolatorTypeIn),
110 electricArrayReflectionTypeSet(BDSArrayReflectionTypeSet()),
111 cacheTransforms(cacheTransformsIn),
112 eScaling(eScalingIn),
113 bScaling(bScalingIn),
114 timeOffset(timeOffsetIn),
115 autoScale(autoScaleIn),
116 stepLimit(stepLimitIn),
117 poleTipRadius(poleTipRadiusIn),
118 beamPipeRadius(beamPipeRadiusIn),
119 chordStepMinimum(-1),
120 secondFieldOnLeft(leftIn),
121 magneticSubFieldName(magneticSubFieldNameIn),
122 electricSubFieldName(electricSubFieldNameIn),
123 usePlacementWorldTransform(false),
124 modulatorInfo(nullptr),
125 transformBeamline(nullptr),
126 nameOfParserDefinition("")
127{
128 if (transformIn != G4Transform3D::Identity)
129 {transform = new G4Transform3D(transformIn);}
130
131 // back calculate tilt angle from field transform
132 G4ThreeVector unitY(0,1,0);
133 G4ThreeVector unitYR = unitY.transform(transformIn.getRotation());
134 tilt = -(CLHEP::halfpi - unitYR.getPhi()); // halfpi offset for which axis we choose by convention
135}
136
137BDSFieldInfo::~BDSFieldInfo()
138{
139 delete magnetStrength;
140 delete transform;
141 delete stepLimit;
142 delete transformBeamline;
143}
144
146 fieldType(other.fieldType),
147 brho(other.brho),
148 integratorType(other.integratorType),
149 provideGlobalTransform(other.provideGlobalTransform),
150 transform(nullptr),
151 magneticFieldFilePath(other.magneticFieldFilePath),
152 magneticFieldFormat(other.magneticFieldFormat),
153 magneticInterpolatorType(other.magneticInterpolatorType),
154 magneticArrayReflectionTypeSet(other.magneticArrayReflectionTypeSet),
155 electricFieldFilePath(other.electricFieldFilePath),
156 electricFieldFormat(other.electricFieldFormat),
157 electricInterpolatorType(other.electricInterpolatorType),
158 electricArrayReflectionTypeSet(other.electricArrayReflectionTypeSet),
159 cacheTransforms(other.cacheTransforms),
160 eScaling(other.eScaling),
161 bScaling(other.bScaling),
162 timeOffset(other.timeOffset),
163 autoScale(other.autoScale),
164 poleTipRadius(other.poleTipRadius),
165 beamPipeRadius(other.beamPipeRadius),
166 chordStepMinimum(other.chordStepMinimum),
167 tilt(other.tilt),
168 secondFieldOnLeft(other.secondFieldOnLeft),
169 magneticSubFieldName(other.magneticSubFieldName),
170 electricSubFieldName(other.electricSubFieldName),
171 usePlacementWorldTransform(other.usePlacementWorldTransform),
172 modulatorInfo(other.modulatorInfo),
173 transformBeamline(nullptr),
174 nameOfParserDefinition(other.nameOfParserDefinition)
175{
176 if (other.transform)
177 {transform = new G4Transform3D(*other.transform);}
178
179 if (other.magnetStrength)
180 {magnetStrength = new BDSMagnetStrength(*other.magnetStrength);}
181 else
182 {magnetStrength = nullptr;} // also nullptr
183
184 if (other.stepLimit)
185 {stepLimit = new G4UserLimits(*other.stepLimit);}
186 else
187 {stepLimit = nullptr;}
188
189 if (other.transformBeamline)
190 {transformBeamline = new G4Transform3D(*other.transformBeamline);}
191}
192
193void BDSFieldInfo::SetUserLimits(G4UserLimits* userLimitsIn)
194{
195 if (stepLimit != defaultUL)
196 {delete stepLimit;} // shouldn't delete global default step limit!
197 stepLimit = userLimitsIn;
198}
199
201 G4bool warn) const
202{
203 if (stepLimit && (stepLimit != defaultUL))
204 {
205 G4UserLimits* old = stepLimit;
206 stepLimit = BDS::CreateUserLimits(stepLimit, maximumStepSize, 1.0);
207 if (stepLimit == old)
208 {warn = false;} // no change and warning would print out wrong number
209 if ((stepLimit != old) && (old != defaultUL))
210 {delete old;}
211 }
212 else
213 {stepLimit = new G4UserLimits(maximumStepSize);}
214
215 if (warn)
216 {
217 G4cout << "Reducing maximum step size of field definition \"" << nameOfParserDefinition
218 << "\" to " << maximumStepSize << " mm " << G4endl;
219 }
220}
221
222std::ostream& operator<< (std::ostream& out, BDSFieldInfo const& info)
223{
224 out << "Parser definition name: \"" << info.nameOfParserDefinition << "\"" << G4endl;
225 out << "Field type: " << info.fieldType << G4endl;
226 out << "Rigidity: " << info.brho << G4endl;
227 out << "Integrator: " << info.integratorType << G4endl;
228 out << "Global transform? " << info.provideGlobalTransform << G4endl;
229 out << "B map file: " << info.magneticFieldFilePath << G4endl;
230 out << "B map file format: " << info.magneticFieldFormat << G4endl;
231 out << "B interpolator " << info.magneticInterpolatorType << G4endl;
232 out << "B array reflection: " << info.magneticArrayReflectionTypeSet << G4endl;
233 out << "E map file: " << info.electricFieldFilePath << G4endl;
234 out << "E map file format: " << info.electricFieldFormat << G4endl;
235 out << "E interpolator " << info.electricInterpolatorType << G4endl;
236 out << "E array reflection: " << info.electricArrayReflectionTypeSet << G4endl;
237 out << "Transform caching: " << info.cacheTransforms << G4endl;
238 out << "E Scaling: " << info.eScaling << G4endl;
239 out << "B Scaling: " << info.bScaling << G4endl;
240 out << "t offset " << info.timeOffset << G4endl;
241 out << "Auto scale " << info.autoScale << G4endl;
242 out << "Pole tip radius: " << info.poleTipRadius << G4endl;
243 out << "Beam pipe radius: " << info.beamPipeRadius << G4endl;
244 out << "Chord Step Min: " << info.chordStepMinimum << G4endl;
245 out << "Tilt: " << info.tilt << G4endl;
246 out << "Second field on left " << info.secondFieldOnLeft << G4endl;
247 out << "Magnetic Sub Field " << info.magneticSubFieldName << G4endl;
248 out << "Electric Sub Field " << info.electricSubFieldName << G4endl;
249 if (info.modulatorInfo)
250 {out << "Modulator " << *(info.modulatorInfo) << G4endl;}
251 out << "Use Placement World Transform " << info.usePlacementWorldTransform << G4endl;
252
253 if (info.magnetStrength)
254 {out << "Magnet strength: " << *(info.magnetStrength) << G4endl;}
255 if (info.stepLimit)
256 {
257 G4Track t = G4Track(); // dummy track
258 G4double maxStep = info.stepLimit->GetMaxAllowedStep(t);
259 out << "Step limit: " << maxStep << G4endl;
260 }
261 return out;
262}
263
264void BDSFieldInfo::Translate(const G4ThreeVector& translationIn)
265{
266 if (!transform)
267 {transform = new G4Transform3D();}
268 G4RotationMatrix rm = transform->getRotation();
269 G4ThreeVector translation = transform->getTranslation();
270 translation += translationIn;
271 G4Transform3D* newTransform = new G4Transform3D(rm, translation);
272 delete transform;
273 transform = newTransform;
274}
275
276G4Transform3D BDSFieldInfo::Transform() const
277{
278 return transform ? *transform : G4Transform3D::Identity;
279}
280
282{
283 return transformBeamline ? *transformBeamline : G4Transform3D::Identity;
284}
285
287{
288 return Transform() * TransformBeamline();
289}
290
291void BDSFieldInfo::SetTransform(const G4Transform3D& transformIn)
292{
293 delete transform;
294 transform = new G4Transform3D(transformIn);
295}
296
297void BDSFieldInfo::SetTransformBeamline(const G4Transform3D& transformIn)
298{
299 delete transformBeamline;
300 transformBeamline = new G4Transform3D(transformIn);
301}
All info required to build complete field of any type.
Definition: BDSFieldInfo.hh:66
G4double beamPipeRadius
Optional radius of beam pipe.
void UpdateUserLimitsLengthMaximumStepSize(G4double maximumStepSize, G4bool warn=false) const
static G4UserLimits * defaultUL
Cache of default user limits.
G4Transform3D * transformBeamline
Transform from curvilinear frame to this field - ie beam line bit only.
G4Transform3D TransformComplete() const
Compound transform of field + beam line transform.
void SetTransformBeamline(const G4Transform3D &transformIn)
Set the beam line transform.
G4double poleTipRadius
Radius at which point the field will be scaled to.
G4double tilt
Cache of tilt of field.
void SetTransform(const G4Transform3D &transformIn)
Set the field definition transform.
G4Transform3D Transform() const
Transform for the field definition only.
G4bool secondFieldOnLeft
Flag for case of two-beam field - if not left, it's right.
void Translate(const G4ThreeVector &translationIn)
G4Transform3D TransformBeamline() const
Transform from the curvilinear coordinates to the beam line component.
void SetUserLimits(G4UserLimits *userLimitsIn)
Delete and replace the user limits which this class owns (only if not default ul).
G4Transform3D * transform
Transform w.r.t. solid field will be attached to.
BDSFieldInfo()
Default constructor for zero field effectively.
Definition: BDSFieldInfo.cc:40
Efficient storage of magnet strengths.
G4UserLimits * CreateUserLimits(G4UserLimits *defaultUL, G4double length, G4double fraction=1.6)