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