| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25
26 /// \file G4FieldSetup.cxx
27 /// \brief Implementation of the G4FieldSetup
28 ///
29 /// This code was initially developed in Geant
30 /// (https://github.com/vmc-project)
31 /// and adapted to Geant4.
32 ///
33 /// \author I. Hrivnacova; IJCLab, Orsay
34
35 #include "G4FieldSetup.hh"
36 #include "G4FieldSetupMessenger.hh"
37
38 #include "G4Exception.hh"
39 #include "G4LogicalVolume.hh"
40 #include "G4SystemOfUnits.hh"
41
42 #include "G4BogackiShampine23.hh"
43 #include "G4BogackiShampine45.hh"
44 #include "G4CachedMagneticField.hh"
45 #include "G4CashKarpRKF45.hh"
46 #include "G4ChordFinder.hh"
47 #include "G4ClassicalRK4.hh"
48 #include "G4ConstRK4.hh"
49 #include "G4DormandPrince745.hh"
50 #include "G4DormandPrinceRK56.hh"
51 #include "G4DormandPrinceRK78.hh"
52 #include "G4ElectroMagneticField.hh"
53 #include "G4EqEMFieldWithEDM.hh"
54 #include "G4EqEMFieldWithSpin.hh"
55 #include "G4EqMagElectricField.hh"
56 #include "G4ExactHelixStepper.hh"
57 #include "G4ExplicitEuler.hh"
58 #include "G4FSALIntegrationDriver.hh"
59 #include "G4FieldManager.hh"
60 #include "G4HelixExplicitEuler.hh"
61 #include "G4HelixHeum.hh"
62 #include "G4HelixImplicitEuler.hh"
63 #include "G4HelixMixedStepper.hh"
64 #include "G4HelixSimpleRunge.hh"
65 #include "G4ImplicitEuler.hh"
66 #include "G4MagneticField.hh"
67 #include "G4MagErrorStepper.hh"
68 #include "G4MagHelicalStepper.hh"
69 #include "G4MagIntegratorDriver.hh"
70 #include "G4Mag_EqRhs.hh"
71 #include "G4Mag_SpinEqRhs.hh"
72 #include "G4Mag_UsualEqRhs.hh"
73 #include "G4NystromRK4.hh"
74 #include "G4RK547FEq1.hh"
75 #include "G4RK547FEq2.hh"
76 #include "G4RK547FEq3.hh"
77 #include "G4RKG3_Stepper.hh"
78 #include "G4SimpleHeum.hh"
79 #include "G4SimpleRunge.hh"
80 #include "G4TsitourasRK45.hh"
81 #include "G4VIntegrationDriver.hh"
82
83
84 //____________________________________________
85 G4FieldSetup::G4FieldSetup(const G4FieldParame
86 G4Field* field, G4LogicalVolume* lv)
87 : fParameters(parameters), fG4Field(field),
88 {
89 // Standard constructor
90
91 fMessenger = new G4FieldSetupMessenger(this)
92
93 // Get or create field manager
94 if (fLogicalVolume == nullptr) {
95 // global field
96 fFieldManager = G4FieldManager::GetGlobalF
97 }
98 else {
99 // local field
100 fFieldManager = new G4FieldManager();
101 G4bool overwriteDaughtersField = true;
102 // TO DO: this parameter should be made
103 fLogicalVolume->SetFieldManager(fFieldMana
104 }
105 }
106
107 //____________________________________________
108 G4FieldSetup::~G4FieldSetup()
109 {
110 // Destructor
111 delete fG4Field;
112 delete fChordFinder;
113 delete fStepper;
114 }
115
116 //
117 // private methods
118 //
119
120 //____________________________________________
121 G4Field* G4FieldSetup::CreateCachedField(
122 const G4FieldParameters& parameters, G4Fie
123 {
124 // Create cached magnetic field if const dista
125 // and field is of G4MagneticField.
126 // Return the input field otherwise.
127
128 if (parameters.GetConstDistance() > 0.) {
129 auto magField = dynamic_cast<G4MagneticFie
130 if (magField == nullptr) {
131 G4Exception(
132 "G4FieldSetup::CreateCachedField:", "G
133 JustWarning, "Incompatible field type.
134 return field;
135 }
136 return new G4CachedMagneticField(magField,
137 }
138
139 return field;
140 }
141
142 //____________________________________________
143 G4EquationOfMotion* G4FieldSetup::CreateEquati
144 {
145 // Set the equation of motion of a particle
146
147 // magnetic fields
148 G4MagneticField* magField = nullptr;
149 if (equation == kEqMagnetic || equation == k
150 magField = dynamic_cast<G4MagneticField*>(
151 if (magField == nullptr) {
152 G4Exception(
153 "G4FieldSetup::CreateEquation:", "Geom
154 FatalErrorInArgument, "Incompatible fi
155 "The field type must be set explicitly
156 return nullptr;
157 }
158 }
159
160 // electromagnetic fields
161 G4ElectroMagneticField* elMagField = nullptr
162 if (equation >= kEqElectroMagnetic && equati
163 elMagField = dynamic_cast<G4ElectroMagneti
164 if (elMagField == nullptr) {
165 G4Exception(
166 "G4FieldSetup::CreateEquation:", "Geom
167 FatalErrorInArgument, "Incompatible fi
168 "The field type must be set explicitly
169 return nullptr;
170 }
171 }
172
173 // electromagnetic fields
174 switch (equation) {
175 case kEqMagnetic:
176 return new G4Mag_UsualEqRhs(magField);
177 break;
178
179 case kEqMagneticWithSpin:
180 return new G4Mag_SpinEqRhs(magField);
181 break;
182
183 case kEqElectroMagnetic:
184 return new G4EqMagElectricField(elMagFie
185 break;
186
187 case kEqEMfieldWithSpin:
188 return new G4EqEMFieldWithSpin(elMagFiel
189 break;
190
191 case kEqEMfieldWithEDM:
192 return new G4EqEMFieldWithEDM(elMagField
193 break;
194 case kUserEquation:
195 // nothing to be done
196 return nullptr;
197 break;
198 }
199
200 G4Exception(
201 "G4FieldSetup::CreateEquation:", "GeomFiel
202 FatalErrorInArgument, "Unknown equation ty
203 return nullptr;
204 }
205
206 //____________________________________________
207 G4MagIntegratorStepper* G4FieldSetup::CreateSt
208 G4EquationOfMotion* equation, G4StepperType
209 {
210 // Set the integrator of particle's equation
211
212 // Check steppers which require equation of
213 auto eqRhs = dynamic_cast<G4Mag_EqRhs*>(equa
214 if ((eqRhs == nullptr) && (stepper > kTsitou
215 G4Exception(
216 "G4FieldSetup::CreateStepper:", "GeomFie
217 FatalErrorInArgument,
218 "The stepper type requires equation of m
219 return nullptr;
220 }
221
222 switch (stepper) {
223 case kBogackiShampine23:
224 return new G4BogackiShampine23(equation)
225 break;
226
227 case kBogackiShampine45:
228 return new G4BogackiShampine45(equation)
229 break;
230
231 case kCashKarpRKF45:
232 return new G4CashKarpRKF45(equation);
233 break;
234
235 case kClassicalRK4:
236 return new G4ClassicalRK4(equation);
237 break;
238
239 case kDormandPrince745:
240 return new G4DormandPrince745(equation);
241 break;
242
243 case kDormandPrinceRK56:
244 return new G4DormandPrinceRK56(equation)
245 break;
246
247 case kDormandPrinceRK78:
248 return new G4DormandPrinceRK78(equation)
249 break;
250
251 case kExplicitEuler:
252 return new G4ExplicitEuler(equation);
253 break;
254
255 case kImplicitEuler:
256 return new G4ImplicitEuler(equation);
257 break;
258
259 case kSimpleHeum:
260 return new G4SimpleHeum(equation);
261 break;
262
263 case kSimpleRunge:
264 return new G4SimpleRunge(equation);
265 break;
266
267 case kTsitourasRK45:
268 return new G4TsitourasRK45(equation);
269 break;
270
271 case kConstRK4:
272 return new G4ConstRK4(eqRhs);
273 break;
274
275 case kExactHelixStepper:
276 return new G4ExactHelixStepper(eqRhs);
277 break;
278
279 case kHelixExplicitEuler:
280 return new G4HelixExplicitEuler(eqRhs);
281 break;
282
283 case kHelixHeum:
284 return new G4HelixHeum(eqRhs);
285 break;
286
287 case kHelixImplicitEuler:
288 return new G4HelixImplicitEuler(eqRhs);
289 break;
290
291 case kHelixMixedStepper:
292 return new G4HelixMixedStepper(eqRhs);
293 break;
294
295 case kHelixSimpleRunge:
296 return new G4HelixSimpleRunge(eqRhs);
297 break;
298
299 case kNystromRK4:
300 return new G4NystromRK4(eqRhs);
301 break;
302
303 case kRKG3Stepper:
304 return new G4RKG3_Stepper(eqRhs);
305 break;
306 case kUserStepper:
307 // nothing to be done
308 return nullptr;
309 break;
310 default:
311 G4Exception(
312 "G4FieldSetup::CreateStepper:", "GeomF
313 FatalErrorInArgument, "Incorrect stepp
314 return nullptr;
315 }
316 }
317
318 //____________________________________________
319 G4VIntegrationDriver* G4FieldSetup::CreateFSAL
320 G4EquationOfMotion* equation, G4StepperType
321 {
322 // Set the FSAL integrator of particle's equ
323
324 switch (stepperType) {
325 case kRK547FEq1:
326 return new G4FSALIntegrationDriver<G4RK5
327 minStep, new G4RK547FEq1(equation));
328
329 case kRK547FEq2:
330 return new G4FSALIntegrationDriver<G4RK5
331 minStep, new G4RK547FEq2(equation));
332
333 case kRK547FEq3:
334 return new G4FSALIntegrationDriver<G4RK5
335 minStep, new G4RK547FEq3(equation));
336
337 default:
338 G4Exception(
339 "G4FieldSetup::CreateFSALStepperAndDri
340 FatalErrorInArgument, "Incorrect stepp
341 return nullptr;
342 }
343 }
344
345 //____________________________________________
346 void G4FieldSetup::CreateCachedField()
347 {
348 // Create cached field (if ConstDistance is
349 fG4Field = CreateCachedField(fParameters, fG
350 }
351
352 //____________________________________________
353 void G4FieldSetup::CreateStepper()
354 {
355 // Create equation of motion (or get the use
356 if (fParameters.GetEquationType() == kUserEq
357 fEquation = fParameters.GetUserEquationOfM
358 // G4cout << "user equation: " << fEquatio
359 }
360 else {
361 delete fEquation;
362 fEquation = nullptr;
363 fEquation = CreateEquation(fParameters.Get
364 // G4cout << "created equation: " << fEqua
365 }
366 fEquation->SetFieldObj(fG4Field);
367
368 // Create stepper (or get the user one if d
369 if (fParameters.GetStepperType() == kUserSte
370 // User stepper
371 fStepper = fParameters.GetUserStepper();
372 // G4cout << "user stepper: " << fStepper
373 }
374 else if (fParameters.GetStepperType() >= kRK
375 // FSAL stepper
376 delete fDriver;
377 delete fStepper;
378 fDriver = nullptr;
379 fStepper = nullptr;
380 fDriver = CreateFSALStepperAndDriver(
381 fEquation, fParameters.GetStepperType(),
382 // G4cout << "FSAL driver: " << fDriver <<
383 if (fDriver) {
384 fStepper = fDriver->GetStepper();
385 // G4cout << "FSAL stepper: " << fStepper
386 }
387 }
388 else {
389 // Normal stepper
390 delete fStepper;
391 fStepper = nullptr;
392 fStepper = CreateStepper(fEquation, fParam
393 // G4cout << "stepper: " << fStepper << G4
394 }
395 }
396
397 //____________________________________________
398 void G4FieldSetup::CreateChordFinder()
399 {
400 // Chord finder
401 if (fParameters.GetFieldType() == kMagnetic)
402 if (fDriver) {
403 fChordFinder = new G4ChordFinder(fDriver
404 // G4cout << "chord finder from existing
405 }
406 else {
407 // Chord finder
408 fChordFinder = new G4ChordFinder(static_
409 fParameters.GetMinimumStep(), fStepper
410 // G4cout << "chord finder from paramete
411 // << " minStep: " << fParameters.GetM
412 }
413 fChordFinder->SetDeltaChord(fParameters.Ge
414 // G4cout << "set delta chord: " << fParam
415 }
416 else if (fParameters.GetFieldType() == kElec
417 G4MagInt_Driver* intDriver = new G4MagInt_
418 fParameters.GetMinimumStep(), fStepper,
419 // G4cout << "kElectroMagnetic: intDriver
420 if (intDriver) {
421 // Chord finder
422 fChordFinder = new G4ChordFinder(intDriv
423 // G4cout << "chord finder from intDrive
424 }
425 }
426 }
427
428 //____________________________________________
429 void G4FieldSetup::UpdateFieldManager()
430 {
431 fFieldManager->SetChordFinder(fChordFinder);
432 fFieldManager->SetDetectorField(fG4Field);
433
434 // G4cout << "Go to set other parameters: "
435 // << "minEpsilon: " << fParameters.GetMi
436 // << "maxEpsilon: " << fParameters.GetMa
437 // << "deltaOneStep: " << fParameters.Get
438 // << "deltaIntersection: " << fParameter
439 fFieldManager->SetMinimumEpsilonStep(fParame
440 fFieldManager->SetMaximumEpsilonStep(fParame
441 fFieldManager->SetDeltaOneStep(fParameters.G
442 fFieldManager->SetDeltaIntersection(fParamet
443 }
444
445 //
446 // public methods
447 //
448
449 //____________________________________________
450 void G4FieldSetup::Clear()
451 {
452 // First clean up previous state.
453 delete fChordFinder;
454 fChordFinder = nullptr;
455
456 if (fG4Field == nullptr) {
457 // G4cout << "fG4Field == nullptr: go to c
458 delete fEquation;
459 delete fDriver;
460 delete fStepper;
461 delete fChordFinder;
462 fEquation = nullptr;
463 fDriver = nullptr;
464 fStepper = nullptr;
465 fChordFinder = nullptr;
466 fFieldManager->SetChordFinder(fChordFinder
467 fFieldManager->SetDetectorField(fG4Field);
468 }
469 }
470
471 //____________________________________________
472 void G4FieldSetup::Update()
473 {
474 // Update field with new field parameters
475
476 Clear();
477 if (fG4Field == nullptr) {
478 // No further update needed
479 return;
480 }
481
482 CreateCachedField();
483 CreateStepper();
484 CreateChordFinder();
485 UpdateFieldManager();
486 }
487
488 //____________________________________________
489 void G4FieldSetup::PrintInfo(G4int verboseLeve
490 {
491 if (verboseLevel == 0) return;
492
493 auto fieldType = G4FieldParameters::FieldTyp
494 auto isCachedMagneticField = (fParameters.Ge
495 if (fLogicalVolume == nullptr) {
496 fieldType = "Global";
497 }
498 else {
499 fieldType = "Local (in ";
500 fieldType.append(fLogicalVolume->GetName()
501 fieldType.append(")");
502 }
503 if (isCachedMagneticField) {
504 fieldType.append(" cached");
505 }
506
507 G4cout << fieldType << " field " << about <<
508 G4cout << G4FieldParameters::StepperTypeName
509 fParameters.GetStepperType())
510 << G4endl;
511
512 if (verboseLevel > 1) {
513 fParameters.PrintParameters();
514 }
515 }
516