| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 /// \file field/field04/src/F04GlobalField.cc <<
28 /// \brief Implementation of the F04GlobalFiel <<
29 // 27 //
30 28
31 #include "F04GlobalField.hh" << 29 #include <time.h>
32 30
33 #include "F04FocusSolenoid.hh" << 31 #include "Randomize.hh"
34 #include "F04SimpleSolenoid.hh" << 32 #include "G4TransportationManager.hh"
35 33
36 #include "G4CashKarpRKF45.hh" <<
37 #include "G4ClassicalRK4.hh" <<
38 #include "G4ExplicitEuler.hh" 34 #include "G4ExplicitEuler.hh"
39 #include "G4ImplicitEuler.hh" 35 #include "G4ImplicitEuler.hh"
40 #include "G4SimpleHeum.hh" <<
41 #include "G4SimpleRunge.hh" 36 #include "G4SimpleRunge.hh"
42 #include "G4SystemOfUnits.hh" << 37 #include "G4SimpleHeum.hh"
43 #include "G4TransportationManager.hh" << 38 #include "G4ClassicalRK4.hh"
44 #include "Randomize.hh" << 39 #include "G4CashKarpRKF45.hh"
45 <<
46 #include <ctime> <<
47 <<
48 //....oooOO0OOooo........oooOO0OOooo........oo <<
49 40
50 G4ThreadLocal F04GlobalField* F04GlobalField:: << 41 #include "F04GlobalField.hh"
51 42
52 //....oooOO0OOooo........oooOO0OOooo........oo << 43 F04GlobalField* F04GlobalField::object = 0;
53 44
54 F04GlobalField::F04GlobalField(F04DetectorCons << 45 F04GlobalField::F04GlobalField() : G4ElectroMagneticField(),
>> 46 minStep(0.01*mm), deltaChord(3.0*mm),
>> 47 deltaOneStep(0.01*mm), deltaIntersection(0.1*mm),
>> 48 epsMin(2.5e-7*mm), epsMax(0.05*mm),
>> 49 fEquation(0), fFieldManager(0),
>> 50 fFieldPropagator(0), fStepper(0), fChordFinder(0)
>> 51 //F04GlobalField::F04GlobalField() : G4MagneticField(),
>> 52 // minStep(0.01*mm), deltaChord(3.0*mm),
>> 53 // deltaOneStep(0.01*mm), deltaIntersection(0.1*mm),
>> 54 // epsMin(2.5e-7*mm), epsMax(0.05*mm),
>> 55 // fEquation(0), fFieldManager(0),
>> 56 // fFieldPropagator(0), fStepper(0), fChordFinder(0)
55 { 57 {
56 fFieldMessenger = new F04FieldMessenger(this << 58 fFieldMessenger = new F04FieldMessenger(this);
57 59
58 fFields = new FieldList(); << 60 fields = new FieldList();
>> 61
>> 62 fStepperType = 4 ; // ClassicalRK4 is default stepper
59 63
60 // set object 64 // set object
61 fObject = this; <<
62 65
63 ConstructField(); << 66 object = this;
64 } <<
65 67
66 //....oooOO0OOooo........oooOO0OOooo........oo << 68 updateField();
>> 69 }
67 70
68 F04GlobalField::~F04GlobalField() 71 F04GlobalField::~F04GlobalField()
69 { 72 {
70 Clear(); << 73 clear();
71 <<
72 delete fFields; <<
73 74
74 delete fFieldMessenger; 75 delete fFieldMessenger;
75 76
76 delete fEquation; << 77 if (fEquation) delete fEquation;
77 delete fStepper; << 78 if (fFieldManager) delete fFieldManager;
78 delete fChordFinder; << 79 if (fFieldPropagator) delete fFieldPropagator;
>> 80 if (fStepper) delete fStepper;
>> 81 if (fChordFinder) delete fChordFinder;
79 } 82 }
80 83
81 //....oooOO0OOooo........oooOO0OOooo........oo << 84 void F04GlobalField::updateField()
82 <<
83 void F04GlobalField::ConstructField() <<
84 { 85 {
85 Clear(); << 86 first = true;
>> 87
>> 88 nfp = 0;
>> 89 fp = 0;
>> 90
>> 91 clear();
86 92
87 // Construct equ. of motion of particles th 93 // Construct equ. of motion of particles through B fields
88 // fEquation = new G4Mag_EqRhs(this); << 94 // fEquation = new G4Mag_EqRhs(this);
89 // Construct equ. of motion of particles th 95 // Construct equ. of motion of particles through e.m. fields
90 // fEquation = new G4EqMagElectricField(thi << 96 // fEquation = new G4EqMagElectricField(this);
91 // Construct equ. of motion of particles in 97 // Construct equ. of motion of particles including spin through B fields
92 // fEquation = new G4Mag_SpinEqRhs(this); << 98 // fEquation = new G4Mag_SpinEqRhs(this);
93 // Construct equ. of motion of particles in 99 // Construct equ. of motion of particles including spin through e.m. fields
94 fEquation = new G4EqEMFieldWithSpin(this); 100 fEquation = new G4EqEMFieldWithSpin(this);
95 101
96 // Get transportation, field, and propagato 102 // Get transportation, field, and propagator managers
97 G4TransportationManager* transportManager = << 103 G4TransportationManager* fTransportManager =
>> 104 G4TransportationManager::GetTransportationManager();
98 105
99 fFieldManager = GetGlobalFieldManager(); 106 fFieldManager = GetGlobalFieldManager();
100 107
101 fFieldPropagator = transportManager->GetProp << 108 fFieldPropagator = fTransportManager->GetPropagatorInField();
102 109
103 // Need to SetFieldChangesEnergy to account 110 // Need to SetFieldChangesEnergy to account for a time varying electric
104 // field (r.f. fields) 111 // field (r.f. fields)
105 fFieldManager->SetFieldChangesEnergy(true); 112 fFieldManager->SetFieldChangesEnergy(true);
106 113
107 // Set the field 114 // Set the field
108 fFieldManager->SetDetectorField(this); 115 fFieldManager->SetDetectorField(this);
109 116
110 // Choose a stepper for integration of the 117 // Choose a stepper for integration of the equation of motion
111 SetStepper(); 118 SetStepper();
112 119
113 // Create a cord finder providing the (glob 120 // Create a cord finder providing the (global field, min step length,
114 // a pointer to the stepper) 121 // a pointer to the stepper)
115 fChordFinder = new G4ChordFinder((G4Magnetic << 122 fChordFinder = new G4ChordFinder((G4MagneticField*)this,minStep,fStepper);
116 123
117 // Set accuracy parameters 124 // Set accuracy parameters
118 fChordFinder->SetDeltaChord(fDeltaChord); << 125 fChordFinder->SetDeltaChord( deltaChord );
119 126
120 fFieldManager->SetAccuraciesWithDeltaOneStep << 127 fFieldManager->SetAccuraciesWithDeltaOneStep(deltaOneStep);
121 128
122 fFieldManager->SetDeltaIntersection(fDeltaIn << 129 fFieldManager->SetDeltaIntersection(deltaIntersection);
123 130
124 fFieldPropagator->SetMinimumEpsilonStep(fEps << 131 fFieldPropagator->SetMinimumEpsilonStep(epsMin);
125 fFieldPropagator->SetMaximumEpsilonStep(fEps << 132 fFieldPropagator->SetMaximumEpsilonStep(epsMax);
126 133
127 G4cout << "Accuracy Parameters:" << 134 G4cout << "Accuracy Parameters:" <<
128 << " MinStep=" << fMinStep << " Delta << 135 " MinStep=" << minStep <<
129 << " DeltaOneStep=" << fDeltaOneStep << 136 " DeltaChord=" << deltaChord <<
130 G4cout << " " << 137 " DeltaOneStep=" << deltaOneStep << G4endl;
131 << " DeltaIntersection=" << fDeltaInt << 138 G4cout << " " <<
132 << " EpsMax=" << fEpsMax << G4endl; << 139 " DeltaIntersection=" << deltaIntersection <<
>> 140 " EpsMin=" << epsMin <<
>> 141 " EpsMax=" << epsMax << G4endl;
133 142
134 fFieldManager->SetChordFinder(fChordFinder); 143 fFieldManager->SetChordFinder(fChordFinder);
135 144
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 } 145 }
175 146
176 //....oooOO0OOooo........oooOO0OOooo........oo << 147 F04GlobalField* F04GlobalField::getObject()
177 <<
178 F04GlobalField* F04GlobalField::GetObject() <<
179 { 148 {
180 if (fObject) return fObject; << 149 if (!object) new F04GlobalField();
181 return nullptr; << 150 return object;
182 } 151 }
183 152
184 //....oooOO0OOooo........oooOO0OOooo........oo <<
185 <<
186 void F04GlobalField::SetStepper() 153 void F04GlobalField::SetStepper()
187 { 154 {
188 delete fStepper; << 155 if(fStepper) delete fStepper;
189 156
190 switch (fStepperType) { << 157 switch ( fStepperType )
>> 158 {
191 case 0: 159 case 0:
192 // fStepper = new G4ExplicitEuler( << 160 // fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking
193 fStepper = new G4ExplicitEuler(fEquation << 161 fStepper = new G4ExplicitEuler( fEquation, 12 ); // with spin tracking
194 G4cout << "G4ExplicitEuler is called" << 162 G4cout << "G4ExplicitEuler is called" << G4endl;
195 break; 163 break;
196 case 1: 164 case 1:
197 // fStepper = new G4ImplicitEuler( << 165 // fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking
198 fStepper = new G4ImplicitEuler(fEquation << 166 fStepper = new G4ImplicitEuler( fEquation, 12 ); // with spin tracking
199 G4cout << "G4ImplicitEuler is called" << 167 G4cout << "G4ImplicitEuler is called" << G4endl;
200 break; 168 break;
201 case 2: 169 case 2:
202 // fStepper = new G4SimpleRunge( fE << 170 // fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking
203 fStepper = new G4SimpleRunge(fEquation, << 171 fStepper = new G4SimpleRunge( fEquation, 12 ); // with spin tracking
204 G4cout << "G4SimpleRunge is called" << G 172 G4cout << "G4SimpleRunge is called" << G4endl;
205 break; 173 break;
206 case 3: 174 case 3:
207 // fStepper = new G4SimpleHeum( fEq << 175 // fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking
208 fStepper = new G4SimpleHeum(fEquation, 1 << 176 fStepper = new G4SimpleHeum( fEquation, 12 ); // with spin tracking
209 G4cout << "G4SimpleHeum is called" << G4 177 G4cout << "G4SimpleHeum is called" << G4endl;
210 break; 178 break;
211 case 4: 179 case 4:
212 // fStepper = new G4ClassicalRK4( f << 180 // fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking
213 fStepper = new G4ClassicalRK4(fEquation, << 181 fStepper = new G4ClassicalRK4( fEquation, 12 ); // with spin tracking
214 G4cout << "G4ClassicalRK4 (default) is c 182 G4cout << "G4ClassicalRK4 (default) is called" << G4endl;
215 break; 183 break;
216 case 5: 184 case 5:
217 // fStepper = new G4CashKarpRKF45( << 185 // fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking
218 fStepper = new G4CashKarpRKF45(fEquation << 186 fStepper = new G4CashKarpRKF45( fEquation, 12 ); // with spin tracking
219 G4cout << "G4CashKarpRKF45 is called" << 187 G4cout << "G4CashKarpRKF45 is called" << G4endl;
220 break; 188 break;
221 default: << 189 default: fStepper = 0;
222 fStepper = nullptr; <<
223 } 190 }
224 } 191 }
225 192
226 //....oooOO0OOooo........oooOO0OOooo........oo <<
227 <<
228 G4FieldManager* F04GlobalField::GetGlobalField 193 G4FieldManager* F04GlobalField::GetGlobalFieldManager()
229 { 194 {
230 return G4TransportationManager::GetTransport << 195 return G4TransportationManager::GetTransportationManager()
>> 196 ->GetFieldManager();
231 } 197 }
232 198
233 //....oooOO0OOooo........oooOO0OOooo........oo <<
234 <<
235 void F04GlobalField::GetFieldValue(const G4dou 199 void F04GlobalField::GetFieldValue(const G4double* point, G4double* field) const
236 { 200 {
237 // NOTE: this routine dominates the CPU time 201 // NOTE: this routine dominates the CPU time for tracking.
238 // Using the simple array fFp[] instead of f << 202 // Using the simple array fp[] instead of fields[]
239 // directly sped it up 203 // directly sped it up
240 204
241 field[0] = field[1] = field[2] = field[3] = 205 field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
242 206
243 // protect against Geant4 bug that calls us 207 // protect against Geant4 bug that calls us with point[] NaN.
244 if (point[0] != point[0]) return; << 208 if(point[0] != point[0]) return;
245 209
246 // (can't use fNfp or fFp, as they may chang << 210 // (can't use nfp or fp, as they may change)
247 if (fFirst) ((F04GlobalField*)this)->SetupAr << 211 if (first) ((F04GlobalField*)this)->setupArray(); // (cast away const)
248 212
249 for (int i = 0; i < fNfp; ++i) { << 213 for (int i=0; i<nfp; ++i) {
250 const F04ElementField* p = fFp[i]; << 214 const F04ElementField* p = fp[i];
251 if (p->IsInBoundingBox(point)) { << 215 if (p->isInBoundingBox(point)) {
252 p->AddFieldValue(point, field); << 216 p->addFieldValue(point,field);
253 } << 217 }
254 } 218 }
255 } <<
256 219
257 //....oooOO0OOooo........oooOO0OOooo........oo << 220 }
258 221
259 void F04GlobalField::Clear() << 222 void F04GlobalField::clear()
260 { 223 {
261 if (fFields) { << 224 if (fields) {
262 if (fFields->size() > 0) { << 225 if (fields->size()>0) {
263 FieldList::iterator i; << 226 FieldList::iterator i;
264 for (i = fFields->begin(); i != fFields- << 227 for (i=fields->begin(); i!=fields->end(); ++i) delete *i;
265 delete *i; << 228 fields->clear();
266 fFields->clear(); << 229 }
267 } <<
268 } 230 }
269 231
270 delete[] fFp; << 232 if (fp) delete[] fp;
271 fFirst = true; <<
272 fNfp = 0; <<
273 fFp = nullptr; <<
274 } <<
275 233
276 //....oooOO0OOooo........oooOO0OOooo........oo << 234 first = true;
>> 235
>> 236 nfp = 0;
>> 237 fp = NULL;
>> 238 }
277 239
278 void F04GlobalField::SetupArray() << 240 void F04GlobalField::setupArray()
279 { 241 {
280 fFirst = false; << 242 first = false;
281 fNfp = fFields->size(); << 243 nfp = fields->size();
282 fFp = new const F04ElementField*[fNfp + 1]; << 244 fp = new const F04ElementField* [nfp+1]; // add 1 so it's never 0
283 for (int i = 0; i < fNfp; ++i) << 245 for (int i=0; i<nfp; ++i) fp[i] = (*fields)[i];
284 fFp[i] = (*fFields)[i]; <<
285 } 246 }
286 247