| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 /// \file G4FieldSetup.cxx
27 /// \brief Implementation of the G4FieldSetup class
28 ///
29 /// This code was initially developed in Geant4 VMC package
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 G4FieldParameters& parameters,
86 G4Field* field, G4LogicalVolume* lv)
87 : fParameters(parameters), fG4Field(field), fLogicalVolume(lv)
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::GetGlobalFieldManager();
97 }
98 else {
99 // local field
100 fFieldManager = new G4FieldManager();
101 G4bool overwriteDaughtersField = true;
102 // TO DO: this parameter should be made optional for users
103 fLogicalVolume->SetFieldManager(fFieldManager, overwriteDaughtersField);
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, G4Field* field)
123 {
124 // Create cached magnetic field if const distance is set > 0.
125 // and field is of G4MagneticField.
126 // Return the input field otherwise.
127
128 if (parameters.GetConstDistance() > 0.) {
129 auto magField = dynamic_cast<G4MagneticField*>(field);
130 if (magField == nullptr) {
131 G4Exception(
132 "G4FieldSetup::CreateCachedField:", "GeomFieldParameters0001",
133 JustWarning, "Incompatible field type.");
134 return field;
135 }
136 return new G4CachedMagneticField(magField, parameters.GetConstDistance());
137 }
138
139 return field;
140 }
141
142 //_____________________________________________________________________________
143 G4EquationOfMotion* G4FieldSetup::CreateEquation(G4EquationType equation)
144 {
145 // Set the equation of motion of a particle in a field
146
147 // magnetic fields
148 G4MagneticField* magField = nullptr;
149 if (equation == kEqMagnetic || equation == kEqMagneticWithSpin) {
150 magField = dynamic_cast<G4MagneticField*>(fG4Field);
151 if (magField == nullptr) {
152 G4Exception(
153 "G4FieldSetup::CreateEquation:", "GeomFieldParameters0001",
154 FatalErrorInArgument, "Incompatible field and equation.\n"
155 "The field type must be set explicitly for other than magnetic field.");
156 return nullptr;
157 }
158 }
159
160 // electromagnetic fields
161 G4ElectroMagneticField* elMagField = nullptr;
162 if (equation >= kEqElectroMagnetic && equation <= kEqEMfieldWithEDM) {
163 elMagField = dynamic_cast<G4ElectroMagneticField*>(fG4Field);
164 if (elMagField == nullptr) {
165 G4Exception(
166 "G4FieldSetup::CreateEquation:", "GeomFieldParameters0001",
167 FatalErrorInArgument, "Incompatible field and equation.\n"
168 "The field type must be set explicitly for other than magnetic field.");
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(elMagField);
185 break;
186
187 case kEqEMfieldWithSpin:
188 return new G4EqEMFieldWithSpin(elMagField);
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:", "GeomFieldParameters0001",
202 FatalErrorInArgument, "Unknown equation type.");
203 return nullptr;
204 }
205
206 //_____________________________________________________________________________
207 G4MagIntegratorStepper* G4FieldSetup::CreateStepper(
208 G4EquationOfMotion* equation, G4StepperType stepper)
209 {
210 // Set the integrator of particle's equation of motion
211
212 // Check steppers which require equation of motion of G4Mag_EqRhs type
213 auto eqRhs = dynamic_cast<G4Mag_EqRhs*>(equation);
214 if ((eqRhs == nullptr) && (stepper > kTsitourasRK45)) {
215 G4Exception(
216 "G4FieldSetup::CreateStepper:", "GeomFieldParameters0001",
217 FatalErrorInArgument,
218 "The stepper type requires equation of motion of G4Mag_EqRhs type.");
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:", "GeomFieldParameters0001",
313 FatalErrorInArgument, "Incorrect stepper type.");
314 return nullptr;
315 }
316 }
317
318 //_____________________________________________________________________________
319 G4VIntegrationDriver* G4FieldSetup::CreateFSALStepperAndDriver(
320 G4EquationOfMotion* equation, G4StepperType stepperType, G4double minStep)
321 {
322 // Set the FSAL integrator of particle's equation of motion
323
324 switch (stepperType) {
325 case kRK547FEq1:
326 return new G4FSALIntegrationDriver<G4RK547FEq1>(
327 minStep, new G4RK547FEq1(equation));
328
329 case kRK547FEq2:
330 return new G4FSALIntegrationDriver<G4RK547FEq2>(
331 minStep, new G4RK547FEq2(equation));
332
333 case kRK547FEq3:
334 return new G4FSALIntegrationDriver<G4RK547FEq3>(
335 minStep, new G4RK547FEq3(equation));
336
337 default:
338 G4Exception(
339 "G4FieldSetup::CreateFSALStepperAndDriver", "GeomFieldParameters0001",
340 FatalErrorInArgument, "Incorrect stepper type.");
341 return nullptr;
342 }
343 }
344
345 //_____________________________________________________________________________
346 void G4FieldSetup::CreateCachedField()
347 {
348 // Create cached field (if ConstDistance is set)
349 fG4Field = CreateCachedField(fParameters, fG4Field);
350 }
351
352 //_____________________________________________________________________________
353 void G4FieldSetup::CreateStepper()
354 {
355 // Create equation of motion (or get the user one if defined)
356 if (fParameters.GetEquationType() == kUserEquation) {
357 fEquation = fParameters.GetUserEquationOfMotion();
358 // G4cout << "user equation: " << fEquation << G4endl;
359 }
360 else {
361 delete fEquation;
362 fEquation = nullptr;
363 fEquation = CreateEquation(fParameters.GetEquationType());
364 // G4cout << "created equation: " << fEquation << G4endl;
365 }
366 fEquation->SetFieldObj(fG4Field);
367
368 // Create stepper (or get the user one if defined)
369 if (fParameters.GetStepperType() == kUserStepper) {
370 // User stepper
371 fStepper = fParameters.GetUserStepper();
372 // G4cout << "user stepper: " << fStepper << G4endl;
373 }
374 else if (fParameters.GetStepperType() >= kRK547FEq1) {
375 // FSAL stepper
376 delete fDriver;
377 delete fStepper;
378 fDriver = nullptr;
379 fStepper = nullptr;
380 fDriver = CreateFSALStepperAndDriver(
381 fEquation, fParameters.GetStepperType(), fParameters.GetMinimumStep());
382 // G4cout << "FSAL driver: " << fDriver << G4endl;
383 if (fDriver) {
384 fStepper = fDriver->GetStepper();
385 // G4cout << "FSAL stepper: " << fStepper << G4endl;
386 }
387 }
388 else {
389 // Normal stepper
390 delete fStepper;
391 fStepper = nullptr;
392 fStepper = CreateStepper(fEquation, fParameters.GetStepperType());
393 // G4cout << "stepper: " << fStepper << G4endl;
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 driver: " << fDriver << G4endl;
405 }
406 else {
407 // Chord finder
408 fChordFinder = new G4ChordFinder(static_cast<G4MagneticField*>(fG4Field),
409 fParameters.GetMinimumStep(), fStepper);
410 // G4cout << "chord finder from parameters: " << fChordFinder
411 // << " minStep: " << fParameters.GetMinimumStep() << G4endl;
412 }
413 fChordFinder->SetDeltaChord(fParameters.GetDeltaChord());
414 // G4cout << "set delta chord: " << fParameters.GetDeltaChord() << G4endl;
415 }
416 else if (fParameters.GetFieldType() == kElectroMagnetic) {
417 G4MagInt_Driver* intDriver = new G4MagInt_Driver(
418 fParameters.GetMinimumStep(), fStepper, fStepper->GetNumberOfVariables());
419 // G4cout << "kElectroMagnetic: intDriver " << intDriver << G4endl;
420 if (intDriver) {
421 // Chord finder
422 fChordFinder = new G4ChordFinder(intDriver);
423 // G4cout << "chord finder from intDriver " << fChordFinder << G4endl;
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.GetMinimumEpsilonStep() << ", "
436 // << "maxEpsilon: " << fParameters.GetMaximumEpsilonStep() << ", "
437 // << "deltaOneStep: " << fParameters.GetDeltaOneStep() << ", "
438 // << "deltaIntersection: " << fParameters.GetDeltaIntersection() << G4endl;
439 fFieldManager->SetMinimumEpsilonStep(fParameters.GetMinimumEpsilonStep());
440 fFieldManager->SetMaximumEpsilonStep(fParameters.GetMaximumEpsilonStep());
441 fFieldManager->SetDeltaOneStep(fParameters.GetDeltaOneStep());
442 fFieldManager->SetDeltaIntersection(fParameters.GetDeltaIntersection());
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 clean-up " << G4endl;
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 verboseLevel, const G4String about)
490 {
491 if (verboseLevel == 0) return;
492
493 auto fieldType = G4FieldParameters::FieldTypeName(fParameters.GetFieldType());
494 auto isCachedMagneticField = (fParameters.GetConstDistance() > 0.);
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 << " with stepper ";
508 G4cout << G4FieldParameters::StepperTypeName(
509 fParameters.GetStepperType())
510 << G4endl;
511
512 if (verboseLevel > 1) {
513 fParameters.PrintParameters();
514 }
515 }
516