Geant4  10.03
F04GlobalField.cc
Go to the documentation of this file.
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 // $Id: F04GlobalField.cc 96981 2016-05-19 09:42:40Z gcosmo $
27 //
30 //
31 
32 #include <time.h>
33 
34 #include "Randomize.hh"
36 
37 #include "G4ExplicitEuler.hh"
38 #include "G4ImplicitEuler.hh"
39 #include "G4SimpleRunge.hh"
40 #include "G4SimpleHeum.hh"
41 #include "G4ClassicalRK4.hh"
42 #include "G4CashKarpRKF45.hh"
43 #include "G4SystemOfUnits.hh"
44 
45 #include "F04GlobalField.hh"
46 #include "F04SimpleSolenoid.hh"
47 #include "F04FocusSolenoid.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 
57  fMinStep(0.01*mm), fDeltaChord(3.0*mm),
58  fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm),
59  fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm),
60  fEquation(0), fFieldManager(0),
61  fFieldPropagator(0), fStepper(0), fChordFinder(0),
62  fDetectorConstruction(det)
63 //F04GlobalField::F04GlobalField(F04DetectorConstruction* det)
64 // : G4MagneticField(),
65 // fMinStep(0.01*mm), fDeltaChord(3.0*mm),
66 // fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm),
67 // fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm),
68 // fEquation(0), fFieldManager(0),
69 // fFieldPropagator(0), fStepper(0), fChordFinder(0),
70 // fDetectorConstruction(det)
71 {
72  fFieldMessenger = new F04FieldMessenger(this,det);
73 
74  fFields = new FieldList();
75 
76  fStepperType = 4 ; // ClassicalRK4 is default stepper
77 
78  // set object
79 
80  fObject = this;
81  fFirst = true;
82 
83  fNfp = 0;
84  fFp = NULL;
85 
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
92 {
93  Clear();
94 
95  delete fFields;
96 
97  delete fFieldMessenger;
98 
99  if (fEquation) delete fEquation;
100  if (fStepper) delete fStepper;
101  if (fChordFinder) delete fChordFinder;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107 {
108  Clear();
109 
110  // Construct equ. of motion of particles through B fields
111 // fEquation = new G4Mag_EqRhs(this);
112  // Construct equ. of motion of particles through e.m. fields
113 // fEquation = new G4EqMagElectricField(this);
114  // Construct equ. of motion of particles including spin through B fields
115 // fEquation = new G4Mag_SpinEqRhs(this);
116  // Construct equ. of motion of particles including spin through e.m. fields
117  fEquation = new G4EqEMFieldWithSpin(this);
118 
119  // Get transportation, field, and propagator managers
120  G4TransportationManager* transportManager =
122 
124 
125  fFieldPropagator = transportManager->GetPropagatorInField();
126 
127  // Need to SetFieldChangesEnergy to account for a time varying electric
128  // field (r.f. fields)
130 
131  // Set the field
133 
134  // Choose a stepper for integration of the equation of motion
135  SetStepper();
136 
137  // Create a cord finder providing the (global field, min step length,
138  // a pointer to the stepper)
140 
141  // Set accuracy parameters
143 
145 
147 
150 
151  G4cout << "Accuracy Parameters:" <<
152  " MinStep=" << fMinStep <<
153  " DeltaChord=" << fDeltaChord <<
154  " DeltaOneStep=" << fDeltaOneStep << G4endl;
155  G4cout << " " <<
156  " DeltaIntersection=" << fDeltaIntersection <<
157  " EpsMin=" << fEpsMin <<
158  " EpsMax=" << fEpsMax << G4endl;
159 
161 
162  G4double l = 0.0;
165 
167  G4ThreeVector captureMgntCenter =
169 
170  F04FocusSolenoid* focusSolenoid =
171  new F04FocusSolenoid(B1, B2, l, logicCaptureMgnt,captureMgntCenter);
172  focusSolenoid -> SetHalf(true);
173 
175 
176  G4LogicalVolume* logicTransferMgnt =
178  G4ThreeVector transferMgntCenter =
180 
181  F04SimpleSolenoid* simpleSolenoid =
182  new F04SimpleSolenoid(B, l, logicTransferMgnt,transferMgntCenter);
183 
184  simpleSolenoid->SetColor("1,0,1");
185  simpleSolenoid->SetColor("0,1,1");
186  simpleSolenoid->SetMaxStep(1.5*mm);
187  simpleSolenoid->SetMaxStep(2.5*mm);
188 
189  if (fFields) {
190  if (fFields->size()>0) {
191  FieldList::iterator i;
192  for (i=fFields->begin(); i!=fFields->end(); ++i){
193  (*i)->Construct();
194  }
195  }
196  }
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 
202 {
203  if (!fObject) new F04GlobalField(det);
204  return fObject;
205 }
206 
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208 
210 {
211  if (fObject) return fObject;
212  return NULL;
213 }
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216 
218 {
219  if(fStepper) delete fStepper;
220 
221  switch ( fStepperType )
222  {
223  case 0:
224 // fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking
225  fStepper = new G4ExplicitEuler( fEquation, 12 ); // with spin tracking
226  G4cout << "G4ExplicitEuler is called" << G4endl;
227  break;
228  case 1:
229 // fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking
230  fStepper = new G4ImplicitEuler( fEquation, 12 ); // with spin tracking
231  G4cout << "G4ImplicitEuler is called" << G4endl;
232  break;
233  case 2:
234 // fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking
235  fStepper = new G4SimpleRunge( fEquation, 12 ); // with spin tracking
236  G4cout << "G4SimpleRunge is called" << G4endl;
237  break;
238  case 3:
239 // fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking
240  fStepper = new G4SimpleHeum( fEquation, 12 ); // with spin tracking
241  G4cout << "G4SimpleHeum is called" << G4endl;
242  break;
243  case 4:
244 // fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking
245  fStepper = new G4ClassicalRK4( fEquation, 12 ); // with spin tracking
246  G4cout << "G4ClassicalRK4 (default) is called" << G4endl;
247  break;
248  case 5:
249 // fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking
250  fStepper = new G4CashKarpRKF45( fEquation, 12 ); // with spin tracking
251  G4cout << "G4CashKarpRKF45 is called" << G4endl;
252  break;
253  default: fStepper = 0;
254  }
255 }
256 
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
258 
260 {
262  ->GetFieldManager();
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 
267 void F04GlobalField::GetFieldValue(const G4double* point, G4double* field) const
268 {
269  // NOTE: this routine dominates the CPU time for tracking.
270  // Using the simple array fFp[] instead of fields[]
271  // directly sped it up
272 
273  field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
274 
275  // protect against Geant4 bug that calls us with point[] NaN.
276  if(point[0] != point[0]) return;
277 
278  // (can't use fNfp or fFp, as they may change)
279  if (fFirst) ((F04GlobalField*)this)->SetupArray(); // (cast away const)
280 
281  for (int i=0; i<fNfp; ++i) {
282  const F04ElementField* p = fFp[i];
283  if (p->IsInBoundingBox(point)) {
284  p->AddFieldValue(point,field);
285  }
286  }
287 
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291 
293 {
294  if (fFields) {
295  if (fFields->size()>0) {
296  FieldList::iterator i;
297  for (i=fFields->begin(); i!=fFields->end(); ++i) delete *i;
298  fFields->clear();
299  }
300  }
301 
302  if (fFp) { delete [] fFp; }
303  fFirst = true;
304  fNfp = 0;
305  fFp = NULL;
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
311 {
312  fFirst = false;
313  fNfp = fFields->size();
314  fFp = new const F04ElementField* [fNfp+1]; // add 1 so it's never 0
315  for (int i=0; i<fNfp; ++i) fFp[i] = (*fFields)[i];
316 }
void SetMaxStep(G4double stp)
SetMaxStep(G4double) sets the max. step size.
G4LogicalVolume * GetTransferMgnt()
G4PropagatorInField * fFieldPropagator
static constexpr double mm
Definition: G4SIunits.hh:115
virtual ~F04GlobalField()
CLHEP::Hep3Vector G4ThreeVector
G4bool SetDetectorField(G4Field *detectorField)
F04DetectorConstruction * fDetectorConstruction
Definition of the F04GlobalField class.
virtual void AddFieldValue(const G4double point[4], G4double field[6]) const =0
AddFieldValue() will add the field value for this element to field[]. Implementations must be sure to...
double B(double temperature)
virtual void GetFieldValue(const G4double *point, G4double *field) const
GetFieldValue() returns the field value at a given point[]. field is really field[6]: Bx...
#define G4ThreadLocal
Definition: tls.hh:89
void SetChordFinder(G4ChordFinder *aChordFinder)
static G4ThreadLocal F04GlobalField * fObject
G4MagIntegratorStepper * fStepper
void SetColor(G4String c)
SetColor(G4String) sets the color.
void SetAccuraciesWithDeltaOneStep(G4double valDeltaOneStep)
static F04GlobalField * GetObject()
G4GLOB_DLL std::ostream G4cout
std::vector< F04ElementField * > FieldList
G4double fDeltaOneStep
G4FieldManager * GetGlobalFieldManager()
Get the global field manager.
Definition of the F04FocusSolenoid class.
G4ChordFinder * fChordFinder
static G4TransportationManager * GetTransportationManager()
F04FieldMessenger * fFieldMessenger
#define fMinStep
G4FieldManager * GetFieldManager() const
G4LogicalVolume * GetCaptureMgnt()
G4EqEMFieldWithSpin * fEquation
bool IsInBoundingBox(const G4double point[4]) const
IsInBoundingBox() returns true if the point is within the global bounding box - global coordinates...
void ConstructField()
constructs all field tracking objects
void SetDeltaIntersection(G4double valueDintersection)
void SetMinimumEpsilonStep(G4double newEpsMin)
void SetFieldChangesEnergy(G4bool value)
Definition of the F04SimpleSolenoid class.
F04GlobalField(F04DetectorConstruction *const)
const F04ElementField ** fFp
#define G4endl
Definition: G4ios.hh:61
FieldList * fFields
G4double fDeltaChord
double G4double
Definition: G4Types.hh:76
void SetMaximumEpsilonStep(G4double newEpsMax)
G4PropagatorInField * GetPropagatorInField() const
void SetStepper()
Set the Stepper.
G4double fDeltaIntersection
void SetDeltaChord(G4double newval)
void Clear()
Clear() removes all ElementField-s from the global object, and destroys them. Used before the geometr...
G4FieldManager * fFieldManager