ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
F04GlobalField.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file F04GlobalField.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 //
26 //
29 //
30 
31 #include <time.h>
32 
33 #include "Randomize.hh"
35 
36 #include "G4ExplicitEuler.hh"
37 #include "G4ImplicitEuler.hh"
38 #include "G4SimpleRunge.hh"
39 #include "G4SimpleHeum.hh"
40 #include "G4ClassicalRK4.hh"
41 #include "G4CashKarpRKF45.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include "F04GlobalField.hh"
45 #include "F04SimpleSolenoid.hh"
46 #include "F04FocusSolenoid.hh"
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
56  fMinStep(0.01*mm), fDeltaChord(3.0*mm),
57  fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm),
58  fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm),
59  fEquation(0), fFieldManager(0),
60  fFieldPropagator(0), fStepper(0), fChordFinder(0),
61  fDetectorConstruction(det)
62 //F04GlobalField::F04GlobalField(F04DetectorConstruction* det)
63 // : G4MagneticField(),
64 // fMinStep(0.01*mm), fDeltaChord(3.0*mm),
65 // fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm),
66 // fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm),
67 // fEquation(0), fFieldManager(0),
68 // fFieldPropagator(0), fStepper(0), fChordFinder(0),
69 // fDetectorConstruction(det)
70 {
71  fFieldMessenger = new F04FieldMessenger(this,det);
72 
73  fFields = new FieldList();
74 
75  fStepperType = 4 ; // ClassicalRK4 is default stepper
76 
77  // set object
78 
79  fObject = this;
80  fFirst = true;
81 
82  fNfp = 0;
83  fFp = NULL;
84 
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
91 {
92  Clear();
93 
94  delete fFields;
95 
96  delete fFieldMessenger;
97 
98  if (fEquation) delete fEquation;
99  if (fStepper) delete fStepper;
100  if (fChordFinder) delete fChordFinder;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
106 {
107  Clear();
108 
109  // Construct equ. of motion of particles through B fields
110 // fEquation = new G4Mag_EqRhs(this);
111  // Construct equ. of motion of particles through e.m. fields
112 // fEquation = new G4EqMagElectricField(this);
113  // Construct equ. of motion of particles including spin through B fields
114 // fEquation = new G4Mag_SpinEqRhs(this);
115  // Construct equ. of motion of particles including spin through e.m. fields
116  fEquation = new G4EqEMFieldWithSpin(this);
117 
118  // Get transportation, field, and propagator managers
119  G4TransportationManager* transportManager =
121 
123 
124  fFieldPropagator = transportManager->GetPropagatorInField();
125 
126  // Need to SetFieldChangesEnergy to account for a time varying electric
127  // field (r.f. fields)
129 
130  // Set the field
132 
133  // Choose a stepper for integration of the equation of motion
134  SetStepper();
135 
136  // Create a cord finder providing the (global field, min step length,
137  // a pointer to the stepper)
139 
140  // Set accuracy parameters
142 
144 
146 
149 
150  G4cout << "Accuracy Parameters:" <<
151  " MinStep=" << fMinStep <<
152  " DeltaChord=" << fDeltaChord <<
153  " DeltaOneStep=" << fDeltaOneStep << G4endl;
154  G4cout << " " <<
155  " DeltaIntersection=" << fDeltaIntersection <<
156  " EpsMin=" << fEpsMin <<
157  " EpsMax=" << fEpsMax << G4endl;
158 
160 
161  G4double l = 0.0;
164 
166  G4ThreeVector captureMgntCenter =
168 
169  F04FocusSolenoid* focusSolenoid =
170  new F04FocusSolenoid(B1, B2, l, logicCaptureMgnt,captureMgntCenter);
171  focusSolenoid -> SetHalf(true);
172 
174 
175  G4LogicalVolume* logicTransferMgnt =
177  G4ThreeVector transferMgntCenter =
179 
180  F04SimpleSolenoid* simpleSolenoid =
181  new F04SimpleSolenoid(B, l, logicTransferMgnt,transferMgntCenter);
182 
183  simpleSolenoid->SetColor("1,0,1");
184  simpleSolenoid->SetColor("0,1,1");
185  simpleSolenoid->SetMaxStep(1.5*mm);
186  simpleSolenoid->SetMaxStep(2.5*mm);
187 
188  if (fFields) {
189  if (fFields->size()>0) {
190  FieldList::iterator i;
191  for (i=fFields->begin(); i!=fFields->end(); ++i){
192  (*i)->Construct();
193  }
194  }
195  }
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
199 
201 {
202  if (!fObject) new F04GlobalField(det);
203  return fObject;
204 }
205 
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
207 
209 {
210  if (fObject) return fObject;
211  return NULL;
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 
217 {
218  if(fStepper) delete fStepper;
219 
220  switch ( fStepperType )
221  {
222  case 0:
223 // fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking
224  fStepper = new G4ExplicitEuler( fEquation, 12 ); // with spin tracking
225  G4cout << "G4ExplicitEuler is called" << G4endl;
226  break;
227  case 1:
228 // fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking
229  fStepper = new G4ImplicitEuler( fEquation, 12 ); // with spin tracking
230  G4cout << "G4ImplicitEuler is called" << G4endl;
231  break;
232  case 2:
233 // fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking
234  fStepper = new G4SimpleRunge( fEquation, 12 ); // with spin tracking
235  G4cout << "G4SimpleRunge is called" << G4endl;
236  break;
237  case 3:
238 // fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking
239  fStepper = new G4SimpleHeum( fEquation, 12 ); // with spin tracking
240  G4cout << "G4SimpleHeum is called" << G4endl;
241  break;
242  case 4:
243 // fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking
244  fStepper = new G4ClassicalRK4( fEquation, 12 ); // with spin tracking
245  G4cout << "G4ClassicalRK4 (default) is called" << G4endl;
246  break;
247  case 5:
248 // fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking
249  fStepper = new G4CashKarpRKF45( fEquation, 12 ); // with spin tracking
250  G4cout << "G4CashKarpRKF45 is called" << G4endl;
251  break;
252  default: fStepper = 0;
253  }
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
257 
259 {
261  ->GetFieldManager();
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265 
267 {
268  // NOTE: this routine dominates the CPU time for tracking.
269  // Using the simple array fFp[] instead of fields[]
270  // directly sped it up
271 
272  field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
273 
274  // protect against Geant4 bug that calls us with point[] NaN.
275  if(point[0] != point[0]) return;
276 
277  // (can't use fNfp or fFp, as they may change)
278  if (fFirst) ((F04GlobalField*)this)->SetupArray(); // (cast away const)
279 
280  for (int i=0; i<fNfp; ++i) {
281  const F04ElementField* p = fFp[i];
282  if (p->IsInBoundingBox(point)) {
283  p->AddFieldValue(point,field);
284  }
285  }
286 
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
290 
292 {
293  if (fFields) {
294  if (fFields->size()>0) {
295  FieldList::iterator i;
296  for (i=fFields->begin(); i!=fFields->end(); ++i) delete *i;
297  fFields->clear();
298  }
299  }
300 
301  if (fFp) { delete [] fFp; }
302  fFirst = true;
303  fNfp = 0;
304  fFp = NULL;
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 
310 {
311  fFirst = false;
312  fNfp = fFields->size();
313  fFp = new const F04ElementField* [fNfp+1]; // add 1 so it's never 0
314  for (int i=0; i<fNfp; ++i) fFp[i] = (*fFields)[i];
315 }