| 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 //
27 /// \file field/field04/src/F04GlobalField.cc
28 /// \brief Implementation of the F04GlobalFiel
29 //
30
31 #include "F04GlobalField.hh"
32
33 #include "F04FocusSolenoid.hh"
34 #include "F04SimpleSolenoid.hh"
35
36 #include "G4CashKarpRKF45.hh"
37 #include "G4ClassicalRK4.hh"
38 #include "G4ExplicitEuler.hh"
39 #include "G4ImplicitEuler.hh"
40 #include "G4SimpleHeum.hh"
41 #include "G4SimpleRunge.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4TransportationManager.hh"
44 #include "Randomize.hh"
45
46 #include <ctime>
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49
50 G4ThreadLocal F04GlobalField* F04GlobalField::
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 F04GlobalField::F04GlobalField(F04DetectorCons
55 {
56 fFieldMessenger = new F04FieldMessenger(this
57
58 fFields = new FieldList();
59
60 // set object
61 fObject = this;
62
63 ConstructField();
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 F04GlobalField::~F04GlobalField()
69 {
70 Clear();
71
72 delete fFields;
73
74 delete fFieldMessenger;
75
76 delete fEquation;
77 delete fStepper;
78 delete fChordFinder;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 void F04GlobalField::ConstructField()
84 {
85 Clear();
86
87 // Construct equ. of motion of particles th
88 // fEquation = new G4Mag_EqRhs(this);
89 // Construct equ. of motion of particles th
90 // fEquation = new G4EqMagElectricField(thi
91 // Construct equ. of motion of particles in
92 // fEquation = new G4Mag_SpinEqRhs(this);
93 // Construct equ. of motion of particles in
94 fEquation = new G4EqEMFieldWithSpin(this);
95
96 // Get transportation, field, and propagato
97 G4TransportationManager* transportManager =
98
99 fFieldManager = GetGlobalFieldManager();
100
101 fFieldPropagator = transportManager->GetProp
102
103 // Need to SetFieldChangesEnergy to account
104 // field (r.f. fields)
105 fFieldManager->SetFieldChangesEnergy(true);
106
107 // Set the field
108 fFieldManager->SetDetectorField(this);
109
110 // Choose a stepper for integration of the
111 SetStepper();
112
113 // Create a cord finder providing the (glob
114 // a pointer to the stepper)
115 fChordFinder = new G4ChordFinder((G4Magnetic
116
117 // Set accuracy parameters
118 fChordFinder->SetDeltaChord(fDeltaChord);
119
120 fFieldManager->SetAccuraciesWithDeltaOneStep
121
122 fFieldManager->SetDeltaIntersection(fDeltaIn
123
124 fFieldPropagator->SetMinimumEpsilonStep(fEps
125 fFieldPropagator->SetMaximumEpsilonStep(fEps
126
127 G4cout << "Accuracy Parameters:"
128 << " MinStep=" << fMinStep << " Delta
129 << " DeltaOneStep=" << fDeltaOneStep
130 G4cout << " "
131 << " DeltaIntersection=" << fDeltaInt
132 << " EpsMax=" << fEpsMax << G4endl;
133
134 fFieldManager->SetChordFinder(fChordFinder);
135
136 G4double l = 0.0;
137 G4double B1 = fDetectorConstruction->GetCapt
138 G4double B2 = fDetectorConstruction->GetCapt
139
140 G4LogicalVolume* logicCaptureMgnt = fDetecto
141 G4ThreeVector captureMgntCenter = fDetectorC
142
143 auto focusSolenoid = new F04FocusSolenoid(B1
144 focusSolenoid->SetHalf(true);
145
146 G4double B = fDetectorConstruction->GetTrans
147
148 G4LogicalVolume* logicTransferMgnt = fDetect
149 G4ThreeVector transferMgntCenter = fDetector
150
151 auto simpleSolenoid = new F04SimpleSolenoid(
152
153 simpleSolenoid->SetColor("1,0,1");
154 simpleSolenoid->SetColor("0,1,1");
155 simpleSolenoid->SetMaxStep(1.5 * mm);
156 simpleSolenoid->SetMaxStep(2.5 * mm);
157
158 if (fFields) {
159 if (fFields->size() > 0) {
160 FieldList::iterator i;
161 for (i = fFields->begin(); i != fFields-
162 (*i)->Construct();
163 }
164 }
165 }
166 }
167
168 //....oooOO0OOooo........oooOO0OOooo........oo
169
170 F04GlobalField* F04GlobalField::GetObject(F04D
171 {
172 if (!fObject) new F04GlobalField(det);
173 return fObject;
174 }
175
176 //....oooOO0OOooo........oooOO0OOooo........oo
177
178 F04GlobalField* F04GlobalField::GetObject()
179 {
180 if (fObject) return fObject;
181 return nullptr;
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186 void F04GlobalField::SetStepper()
187 {
188 delete fStepper;
189
190 switch (fStepperType) {
191 case 0:
192 // fStepper = new G4ExplicitEuler(
193 fStepper = new G4ExplicitEuler(fEquation
194 G4cout << "G4ExplicitEuler is called" <<
195 break;
196 case 1:
197 // fStepper = new G4ImplicitEuler(
198 fStepper = new G4ImplicitEuler(fEquation
199 G4cout << "G4ImplicitEuler is called" <<
200 break;
201 case 2:
202 // fStepper = new G4SimpleRunge( fE
203 fStepper = new G4SimpleRunge(fEquation,
204 G4cout << "G4SimpleRunge is called" << G
205 break;
206 case 3:
207 // fStepper = new G4SimpleHeum( fEq
208 fStepper = new G4SimpleHeum(fEquation, 1
209 G4cout << "G4SimpleHeum is called" << G4
210 break;
211 case 4:
212 // fStepper = new G4ClassicalRK4( f
213 fStepper = new G4ClassicalRK4(fEquation,
214 G4cout << "G4ClassicalRK4 (default) is c
215 break;
216 case 5:
217 // fStepper = new G4CashKarpRKF45(
218 fStepper = new G4CashKarpRKF45(fEquation
219 G4cout << "G4CashKarpRKF45 is called" <<
220 break;
221 default:
222 fStepper = nullptr;
223 }
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oo
227
228 G4FieldManager* F04GlobalField::GetGlobalField
229 {
230 return G4TransportationManager::GetTransport
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oo
234
235 void F04GlobalField::GetFieldValue(const G4dou
236 {
237 // NOTE: this routine dominates the CPU time
238 // Using the simple array fFp[] instead of f
239 // directly sped it up
240
241 field[0] = field[1] = field[2] = field[3] =
242
243 // protect against Geant4 bug that calls us
244 if (point[0] != point[0]) return;
245
246 // (can't use fNfp or fFp, as they may chang
247 if (fFirst) ((F04GlobalField*)this)->SetupAr
248
249 for (int i = 0; i < fNfp; ++i) {
250 const F04ElementField* p = fFp[i];
251 if (p->IsInBoundingBox(point)) {
252 p->AddFieldValue(point, field);
253 }
254 }
255 }
256
257 //....oooOO0OOooo........oooOO0OOooo........oo
258
259 void F04GlobalField::Clear()
260 {
261 if (fFields) {
262 if (fFields->size() > 0) {
263 FieldList::iterator i;
264 for (i = fFields->begin(); i != fFields-
265 delete *i;
266 fFields->clear();
267 }
268 }
269
270 delete[] fFp;
271 fFirst = true;
272 fNfp = 0;
273 fFp = nullptr;
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oo
277
278 void F04GlobalField::SetupArray()
279 {
280 fFirst = false;
281 fNfp = fFields->size();
282 fFp = new const F04ElementField*[fNfp + 1];
283 for (int i = 0; i < fNfp; ++i)
284 fFp[i] = (*fFields)[i];
285 }
286