ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
F01FieldSetup.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file F01FieldSetup.cc
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 //
28 //
29 //
30 //
31 // User Field setup class implementation.
32 //
33 //
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36 
37 #include "F01FieldSetup.hh"
38 #include "F01FieldMessenger.hh"
39 
40 #include "G4MagneticField.hh"
41 #include "G4UniformMagField.hh"
42 #include "G4FieldManager.hh"
43 #include "G4TransportationManager.hh"
44 #include "G4Mag_UsualEqRhs.hh"
45 #include "G4MagIntegratorStepper.hh"
46 #include "G4ChordFinder.hh"
47 
48 #include "G4ExplicitEuler.hh"
49 #include "G4ImplicitEuler.hh"
50 #include "G4SimpleRunge.hh"
51 #include "G4SimpleHeum.hh"
52 #include "G4ClassicalRK4.hh"
53 #include "G4HelixExplicitEuler.hh"
54 #include "G4HelixImplicitEuler.hh"
55 #include "G4HelixSimpleRunge.hh"
56 #include "G4CashKarpRKF45.hh"
57 #include "G4RKG3_Stepper.hh"
58 #include "G4ConstRK4.hh"
59 #include "G4NystromRK4.hh"
60 #include "G4HelixMixedStepper.hh"
61 #include "G4ExactHelixStepper.hh"
62 
63 // Newest steppers - from Release 10.3-beta (June 2013)
64 #include "G4BogackiShampine23.hh"
65 #include "G4BogackiShampine45.hh"
66 #include "G4DormandPrince745.hh"
67 #include "G4DormandPrinceRK56.hh"
68 #include "G4DormandPrinceRK78.hh"
69 #include "G4TsitourasRK45.hh"
70 
71 #include "G4PhysicalConstants.hh"
72 #include "G4SystemOfUnits.hh"
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 enum EStepperNumber { kDormandPrince45 = 17, kBogackiShampine45= 45, kClassicalRK4 = 4,
77  kNystromRK4 = 13 /*soon 40*/,
78  kDormandPrince56 = 56, kBogackiShampine23= 23, kCashKarp = 8,
79  kDormandPrince78 = 78, kTsitouras45 = 145
80 } ;
81 
82 // Constructors:
83 
84 F01FieldSetup::F01FieldSetup(G4ThreeVector fieldVector,
85  G4int stepperNum,
86  G4bool useFSALstepper )
87  : fMagneticField(new G4UniformMagField(fieldVector)),
88  fUseFSALstepper(useFSALstepper),
89  fStepperType(0),
90  fMinStep(0.)
91 {
92  G4cout << " F01FieldSetup: magnetic field set to Uniform( "
93  << fieldVector << " ) " << G4endl;
94 
95  if( stepperNum == -1000 )
96  {
97  fUseFSALstepper = useFSALstepper;
98  if( !useFSALstepper )
99  fStepperType= 17; // Use Dormand Prince (7) 4/5 as default stepper
100  else
101  fStepperType = 101;
102  }
103  else
104  {
105  fUseFSALstepper = ( stepperNum > 0 );
106  if( stepperNum > 0 )
107  fStepperType = stepperNum;
108  else
109  fStepperType = - stepperNum;
110 
111  }
112 
113  InitialiseAll();
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 F01FieldSetup::F01FieldSetup()
119  : fMagneticField(new G4UniformMagField(G4ThreeVector())),
120  fUseFSALstepper(false),
121  fStepperType(17), // Use Dormand Prince (7) 4/5 as default stepper
122  fMinStep(0.)
123 {
124  G4cout << " F01FieldSetup: magnetic field set to Uniform( 0.0, 0, 0 ) "
125  << G4endl;
126  InitialiseAll();
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
131 void F01FieldSetup::InitialiseAll()
132 {
133  fFieldMessenger = new F01FieldMessenger(this);
134 
135  fEquation = new G4Mag_UsualEqRhs(fMagneticField);
136 
137  fMinStep = 1.0*mm; // minimal step of 1 mm is default
138 
139  fFieldManager = G4TransportationManager::GetTransportationManager()
140  ->GetFieldManager();
141 
142  if( fUseFSALstepper ) {
143  CreateFSALStepperAndChordFinder();
144  }
145  else
146  {
147  CreateStepperAndChordFinder();
148  }
149 
150  G4cout << " 4. Updating Field Manager." << G4endl;
151  fFieldManager->SetChordFinder( fChordFinder );
152  fFieldManager->SetDetectorField(fMagneticField );
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 
157 F01FieldSetup::~F01FieldSetup()
158 {
159  delete fMagneticField;
160  delete fChordFinder;
161  delete fStepper;
162  delete fFieldMessenger;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166 
167 void F01FieldSetup::CreateStepperAndChordFinder()
168 {
169  delete fChordFinder;
170  fChordFinder= nullptr;
171 
172  // Update field
173  G4cout << " F01FieldSetup::CreateStepperAndChordFinder() called. " << G4endl
174  << " 1. Creating Stepper." << G4endl;
175 
176  SetStepper();
177  G4cout<<"The minimal step is equal to "<<fMinStep/mm<<" mm"<<G4endl;
178 
179  G4cout << " 2. Creating ChordFinder." << G4endl;
180  fChordFinder = new G4ChordFinder( fMagneticField, fMinStep,fStepper );
181 
182  G4cout << " 3. Updating Field Manager." << G4endl;
183  fFieldManager->SetChordFinder( fChordFinder );
184  fFieldManager->SetDetectorField(fMagneticField );
185 }
186 
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188 
189 void F01FieldSetup::SetStepper()
190 {
191 // Set stepper according to the stepper type
192 
193  if (fStepper) delete fStepper;
194 
195  switch ( fStepperType )
196  {
197  // The new default in G4 and here ( since G4 10.4 Dec 2017 )
198  case 17:
199  case 457:
200  case 745:
201  fStepper = new G4DormandPrince745( fEquation );
202  G4cout<<"G4DormandPrince745 Stepper is chosen"<<G4endl;
203  break;
204 
205  case 0:
206  fStepper = new G4ExplicitEuler( fEquation );
207  G4cout<<"G4ExplicitEuler is chosen."<<G4endl;
208  break;
209  case 1:
210  fStepper = new G4ImplicitEuler( fEquation );
211  G4cout<<"G4ImplicitEuler is chosen"<<G4endl;
212  break;
213  case 2:
214  fStepper = new G4SimpleRunge( fEquation );
215  G4cout<<"G4SimpleRunge is chosen"<<G4endl;
216  break;
217  case 3:
218  fStepper = new G4SimpleHeum( fEquation );
219  G4cout<<"G4SimpleHeum is chosen"<<G4endl;
220  break;
221  case 4:
222  fStepper = new G4ClassicalRK4( fEquation );
223  G4cout<<"G4ClassicalRK4 (default) is chosen"<<G4endl;
224  break;
225  case 5:
226  fStepper = new G4HelixExplicitEuler( fEquation );
227  G4cout<<"G4HelixExplicitEuler is chosen"<<G4endl;
228  break;
229  case 6:
230  fStepper = new G4HelixImplicitEuler( fEquation );
231  G4cout<<"G4HelixImplicitEuler is chosen"<<G4endl;
232  break;
233  case 7:
234  fStepper = new G4HelixSimpleRunge( fEquation );
235  G4cout<<"G4HelixSimpleRunge is chosen"<<G4endl;
236  break;
237  case 8:
238  fStepper = new G4CashKarpRKF45( fEquation );
239  G4cout<<"G4CashKarpRKF45 is chosen"<<G4endl;
240  break;
241  case 9:
242  fStepper = new G4RKG3_Stepper( fEquation );
243  G4cout<<"G4RKG3_Stepper is chosen"<<G4endl;
244  break;
245  case 10:
246  fStepper = new G4ExactHelixStepper( fEquation );
247  G4cout<<"G4ExactHelixStepper is chosen"<<G4endl;
248  break;
249  case 11:
250  fStepper = new G4HelixMixedStepper( fEquation );
251  G4cout<<"G4HelixMixedStepper is chosen"<<G4endl;
252  break;
253  case 12:
254  fStepper = new G4ConstRK4( fEquation );
255  G4cout<<"G4ConstRK4 Stepper is chosen"<<G4endl;
256  break;
257  case 13:
258  case 40:
259  fStepper = new G4NystromRK4( fEquation );
260  G4cout<<" G4NystromRK4 Stepper is chosen"<<G4endl;
261  break;
262  case 14:
263  case 23:
264  fStepper = new G4BogackiShampine23( fEquation );
265  G4cout<<"G4BogackiShampine23 Stepper is chosen"<<G4endl;
266  break;
267 
268  // Other optimised 4/5th order embedded steppers
269  case 15:
270  case 45:
271  fStepper = new G4BogackiShampine45( fEquation );
272  G4cout<<"G4BogackiShampine45 Stepper is chosen"<<G4endl;
273  break;
274 
275  // case 145:
276  case kTsitouras45:
277  fStepper = new G4TsitourasRK45( fEquation );
278  G4cout<<"G4TsitourasRK45 Stepper is chosen"<<G4endl;
279  break;
280 
281  // Higher order embedded steppers - for very smooth fields
282  case 56:
283  fStepper = new G4DormandPrinceRK56( fEquation );
284  G4cout<<"G4DormandPrinceRK56 Stepper is chosen"<<G4endl;
285  break;
286  case 78:
287  fStepper = new G4DormandPrinceRK78( fEquation );
288  G4cout<<"G4DormandPrinceRK78 Stepper is chosen"<<G4endl;
289  break;
290 
291  default:
292  // fStepper = new G4ClassicalRK4( fEquation );
293  // G4cout<<"G4ClassicalRK4 Stepper (default) is chosen"<<G4endl;
294  fStepper = new G4DormandPrince745( fEquation );
295  G4cout<<"G4DormandPrince745 (default) Stepper is chosen"<<G4endl;
296  break;
297  }
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
301 
302 #include "G4VIntegrationDriver.hh"
303 #include "G4FSALIntegrationDriver.hh"
304 #include "G4RK547FEq1.hh"
305 #include "G4RK547FEq2.hh"
306 #include "G4RK547FEq3.hh"
307 
308 G4VIntegrationDriver*
309 F01FieldSetup::CreateFSALStepperAndDriver()
310 {
311  // using FsalStepperType = G4RK547FEq1;
312  const char *methodName= "F01FieldSetup::CreateFSALStepperAndDriver()";
313  if (fStepper) delete fStepper;
314  fStepper = nullptr;
315 
316  G4cout << " F01FieldSetup::CreateFSALStepperAndDriver() called. " << G4endl;
317  G4cout << " 1. Creating Stepper." << G4endl;
318  // auto fsalStepper = new FsalStepperType( fEquation );
319  G4RK547FEq1* stepper1 = nullptr;
320  G4RK547FEq2* stepper2 = nullptr;
321  G4RK547FEq3* stepper3 = nullptr;
322 
323  G4cout << " 2. Creating FSAL Driver." << G4endl;
324  G4VIntegrationDriver* fsalDriver = nullptr;
325  switch ( fStepperType )
326  {
327  case 1:
328  case 101:
329  stepper1 = new G4RK547FEq1( fEquation );
330  fsalDriver = new G4FSALIntegrationDriver<G4RK547FEq1>( fMinStep, stepper1 );
331  G4cout << " Stepper type '1' is G4RK547FEq1 stepper (in FSAL mode) with FSAL driver. "
332  << G4endl;
333  fStepper = stepper1;
334  stepper1 = nullptr;
335  break;
336 
337  case 2:
338  case 102:
339  stepper2= new G4RK547FEq2( fEquation );
340  fsalDriver = new G4FSALIntegrationDriver<G4RK547FEq2>( fMinStep, stepper2 );
341  G4cout << " Stepper type '2' is G4RK547FEq2 stepper (in FSAL mode) with FSAL driver. "
342  << G4endl;
343  fStepper = stepper2;
344  stepper2 = nullptr;
345  break;
346 
347  case 3:
348  case 103:
349  stepper3 = new G4RK547FEq3( fEquation );
350  fsalDriver = new G4FSALIntegrationDriver<G4RK547FEq3>( fMinStep, stepper3 );
351  G4cout << " Stepper type '3' is G4RK547FEq3 stepper (in FSAL mode) with FSAL driver. "
352  << G4endl;
353  fStepper = stepper3;
354  stepper3 = nullptr;
355  break;
356 
357  default:
358  G4cout << " Warning from " << methodName << " : stepperType (= "
359  << fStepperType << " ) is unknown. " << G4endl
360  << " Using value '1' instead - i.e. G4RK547FEq1 stepper. "
361  << G4endl;
362  stepper1 = new G4RK547FEq1( fEquation );
363  fsalDriver = new G4FSALIntegrationDriver<G4RK547FEq1>( fMinStep, stepper1 );
364  fStepper = stepper1;
365  stepper1 = nullptr;
366  break;
367  }
368 
369  delete stepper1; stepper1 = nullptr;
370  delete stepper2; stepper2 = nullptr;
371  delete stepper3; stepper3 = nullptr;
372 
373  if( fsalDriver )
374  fStepper = fsalDriver->GetStepper();
375 
376  return fsalDriver;
377 }
378 
379 void F01FieldSetup::CreateFSALStepperAndChordFinder()
380 {
381  // using FsalStepperType = G4DormandPrince745; // eventually ?
382  delete fChordFinder;
383  fChordFinder= nullptr;
384 
385  G4cout << " F01FieldSetup::CreateFSALStepperAndChordFinder() called. " << G4endl;
386 
387  auto FSALdriver= CreateFSALStepperAndDriver();
388  fDriver = FSALdriver;
389  G4cout<<"The minimal step is equal to "<<fMinStep/mm<<" mm"<<G4endl;
390 
391  G4cout << " 3. Creating ChordFinder." << G4endl;
392  fChordFinder = new G4ChordFinder( FSALdriver ); // ( fMagneticField, fMinStep, fStepper );
393 }
394 
395 
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
397 
398 void F01FieldSetup::SetFieldZValue(G4double fieldStrength)
399 {
400  // Set the value of the Global Field to fieldValue along Z
401 
402  SetFieldValue(G4ThreeVector(0, 0, fieldStrength));
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406 
407 void F01FieldSetup::SetFieldValue(G4ThreeVector fieldVector)
408 {
409  // Set the value of the Global Field
410 
411  if (fMagneticField) delete fMagneticField;
412 
413 #ifdef G4VERBOSE
414  G4cout << "Setting Field strength to "
415  << fieldVector / gauss << " Gauss." << G4endl;
416 #endif
417 
418  if (fieldVector != G4ThreeVector(0.,0.,0.))
419  {
420  fMagneticField = new G4UniformMagField(fieldVector);
421  }
422  else
423  {
424 #ifdef G4VERBOSE
425  G4cout << " Magnetic field pointer is null." << G4endl;
426 #endif
427  // If the new field's value is Zero, signal it as below
428  // so that it is not used for propagation.
429  fMagneticField = 0;
430  }
431 
432  // Set this as the field of the global Field Manager
433  GetGlobalFieldManager()->SetDetectorField(fMagneticField);
434 
435  // Now notify equation of new field
436  fEquation->SetFieldObj( fMagneticField );
437 
438 }
439 
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
441 
442 G4FieldManager* F01FieldSetup::GetGlobalFieldManager()
443 {
444  // Utility method
445 
446  return G4TransportationManager::GetTransportationManager()
447  ->GetFieldManager();
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......