| 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 // G4FieldManager implementation 26 // G4FieldManager implementation
27 // 27 //
28 // Author: John Apostolakis, 10.03.97 - design 28 // Author: John Apostolakis, 10.03.97 - design and implementation
29 // ------------------------------------------- 29 // -------------------------------------------------------------------
30 30
31 #include "G4FieldManager.hh" 31 #include "G4FieldManager.hh"
32 #include "G4Field.hh" 32 #include "G4Field.hh"
33 #include "G4MagneticField.hh" 33 #include "G4MagneticField.hh"
34 #include "G4ChordFinder.hh" 34 #include "G4ChordFinder.hh"
35 #include "G4FieldManagerStore.hh" 35 #include "G4FieldManagerStore.hh"
36 #include "G4SystemOfUnits.hh" 36 #include "G4SystemOfUnits.hh"
37 37
38 G4ThreadLocal G4FieldManager* G4FieldManager:: <<
39 G4double G4FieldManager::fDefault_Delta_One_St 38 G4double G4FieldManager::fDefault_Delta_One_Step_Value= 0.01 * millimeter;
40 G4double G4FieldManager::fDefault_Delta_Inters 39 G4double G4FieldManager::fDefault_Delta_Intersection_Val= 0.001 * millimeter;
41 G4bool G4FieldManager::fVerboseConstruction= <<
42 <<
43 G4double G4FieldManager::fMaxAcceptedEpsilon= <<
44 // Requesting a large epsilon (max) value pr <<
45 // every integration segment. <<
46 // Problems occur because some methods (incl <<
47 // error appears to be a substantial underes <<
48 // So the value for fMaxAcceptedEpsilon is r <<
49 40
50 G4FieldManager::G4FieldManager(G4Field* detect 41 G4FieldManager::G4FieldManager(G4Field* detectorField,
51 G4ChordFinder* 42 G4ChordFinder* pChordFinder,
52 G4bool fieldCha 43 G4bool fieldChangesEnergy )
53 : fDetectorField(detectorField), 44 : fDetectorField(detectorField),
54 fChordFinder(pChordFinder), 45 fChordFinder(pChordFinder),
55 fDelta_One_Step_Value( fDefault_Delta_One 46 fDelta_One_Step_Value( fDefault_Delta_One_Step_Value ),
56 fDelta_Intersection_Val( fDefault_Delta_I 47 fDelta_Intersection_Val( fDefault_Delta_Intersection_Val ),
57 fEpsilonMin( fEpsilonMinDefault ), 48 fEpsilonMin( fEpsilonMinDefault ),
58 fEpsilonMax( fEpsilonMaxDefault) 49 fEpsilonMax( fEpsilonMaxDefault)
59 { 50 {
60 if ( detectorField != nullptr ) 51 if ( detectorField != nullptr )
61 { 52 {
62 fFieldChangesEnergy = detectorField->Does 53 fFieldChangesEnergy = detectorField->DoesFieldChangeEnergy();
63 } 54 }
64 else 55 else
65 { 56 {
66 fFieldChangesEnergy = fieldChangesEnergy; 57 fFieldChangesEnergy = fieldChangesEnergy;
67 } 58 }
68 59
69 if( fVerboseConstruction) <<
70 { <<
71 G4cout << "G4FieldManager/ctor#1 fEpsilon <<
72 } <<
73 <<
74 // Add to store 60 // Add to store
75 // 61 //
76 G4FieldManagerStore::Register(this); 62 G4FieldManagerStore::Register(this);
77 } 63 }
78 64
79 G4FieldManager::G4FieldManager(G4MagneticField 65 G4FieldManager::G4FieldManager(G4MagneticField* detectorField)
80 : fDetectorField(detectorField), fAllocated 66 : fDetectorField(detectorField), fAllocatedChordFinder(true),
81 fDelta_One_Step_Value( fDefault_Delta_O 67 fDelta_One_Step_Value( fDefault_Delta_One_Step_Value ),
82 fDelta_Intersection_Val( fDefault_Delta_I 68 fDelta_Intersection_Val( fDefault_Delta_Intersection_Val ),
83 fEpsilonMin( fEpsilonMinDefault ), 69 fEpsilonMin( fEpsilonMinDefault ),
84 fEpsilonMax( fEpsilonMaxDefault ) 70 fEpsilonMax( fEpsilonMaxDefault )
85 { 71 {
86 fChordFinder = new G4ChordFinder( detectorF 72 fChordFinder = new G4ChordFinder( detectorField );
87 73
88 if( fVerboseConstruction ) <<
89 { <<
90 G4cout << "G4FieldManager/ctor#2 fEpsilon <<
91 } <<
92 // Add to store 74 // Add to store
93 // 75 //
94 G4FieldManagerStore::Register(this); 76 G4FieldManagerStore::Register(this);
95 } 77 }
96 78
97 G4FieldManager* G4FieldManager::Clone() const 79 G4FieldManager* G4FieldManager::Clone() const
98 { 80 {
99 G4Field* aField = nullptr; 81 G4Field* aField = nullptr;
100 G4FieldManager* aFM = nullptr; 82 G4FieldManager* aFM = nullptr;
101 G4ChordFinder* aCF = nullptr; 83 G4ChordFinder* aCF = nullptr;
102 try { 84 try {
103 if ( fDetectorField != nullptr ) 85 if ( fDetectorField != nullptr )
104 { 86 {
105 aField = fDetectorField->Clone(); 87 aField = fDetectorField->Clone();
106 } 88 }
107 89
108 // Create a new field manager, note th 90 // Create a new field manager, note that we do not set
109 // any chordfinder now 91 // any chordfinder now
110 // 92 //
111 aFM = new G4FieldManager( aField , nul 93 aFM = new G4FieldManager( aField , nullptr , fFieldChangesEnergy );
112 94
113 // Check if originally we have the fAl 95 // Check if originally we have the fAllocatedChordFinder variable
114 // set, in case, call chord constructo 96 // set, in case, call chord constructor
115 // 97 //
116 if ( fAllocatedChordFinder ) 98 if ( fAllocatedChordFinder )
117 { 99 {
118 aFM->CreateChordFinder( dynamic_ca 100 aFM->CreateChordFinder( dynamic_cast<G4MagneticField*>(aField) );
119 } 101 }
120 else 102 else
121 { 103 {
122 // Chord was specified by user, sh 104 // Chord was specified by user, should we clone?
123 // TODO: For the moment copy point 105 // TODO: For the moment copy pointer, to be understood
124 // if cloning of ChordFinder 106 // if cloning of ChordFinder is needed
125 // 107 //
126 aCF = fChordFinder; /*->Clone*/ 108 aCF = fChordFinder; /*->Clone*/
127 aFM->fChordFinder = aCF; 109 aFM->fChordFinder = aCF;
128 } 110 }
129 111
130 // Copy values of other variables 112 // Copy values of other variables
131 113
132 aFM->fEpsilonMax = fEpsilonMax; 114 aFM->fEpsilonMax = fEpsilonMax;
133 aFM->fEpsilonMin = fEpsilonMin; 115 aFM->fEpsilonMin = fEpsilonMin;
134 aFM->fDelta_Intersection_Val = fDelta_ 116 aFM->fDelta_Intersection_Val = fDelta_Intersection_Val;
135 aFM->fDelta_One_Step_Value = fDelta_On 117 aFM->fDelta_One_Step_Value = fDelta_One_Step_Value;
136 // TODO: Should we really add to the 118 // TODO: Should we really add to the store the cloned FM?
137 // Who will use this? 119 // Who will use this?
138 } 120 }
139 catch ( ... ) 121 catch ( ... )
140 { 122 {
141 // Failed creating clone: probably use 123 // Failed creating clone: probably user did not implement Clone method
142 // in derived classes? 124 // in derived classes?
143 // Perform clean-up after ourselves... 125 // Perform clean-up after ourselves...
144 delete aField; 126 delete aField;
145 delete aFM; 127 delete aFM;
146 delete aCF; 128 delete aCF;
147 throw; 129 throw;
148 } 130 }
149 <<
150 G4cout << "G4FieldManager/clone fEpsilon M <<
151 return aFM; 131 return aFM;
152 } 132 }
153 133
154 void G4FieldManager::ConfigureForTrack( const 134 void G4FieldManager::ConfigureForTrack( const G4Track * )
155 { 135 {
156 // Default is to do nothing! 136 // Default is to do nothing!
157 } 137 }
158 138
159 G4FieldManager::~G4FieldManager() 139 G4FieldManager::~G4FieldManager()
160 { 140 {
161 if( fAllocatedChordFinder ) 141 if( fAllocatedChordFinder )
162 { 142 {
163 delete fChordFinder; 143 delete fChordFinder;
164 } 144 }
165 G4FieldManagerStore::DeRegister(this); 145 G4FieldManagerStore::DeRegister(this);
166 } 146 }
167 147
168 void 148 void
169 G4FieldManager::CreateChordFinder(G4MagneticFi 149 G4FieldManager::CreateChordFinder(G4MagneticField* detectorMagField)
170 { 150 {
171 if ( fAllocatedChordFinder ) 151 if ( fAllocatedChordFinder )
172 { 152 {
173 delete fChordFinder; 153 delete fChordFinder;
174 } 154 }
175 fAllocatedChordFinder = false; 155 fAllocatedChordFinder = false;
176 156
177 if( detectorMagField != nullptr ) 157 if( detectorMagField != nullptr )
178 { 158 {
179 fChordFinder = new G4ChordFinder( detect 159 fChordFinder = new G4ChordFinder( detectorMagField );
180 fAllocatedChordFinder = true; 160 fAllocatedChordFinder = true;
181 } 161 }
182 else 162 else
183 { 163 {
184 fChordFinder = nullptr; 164 fChordFinder = nullptr;
185 } 165 }
186 } 166 }
187 167
188 void G4FieldManager::InitialiseFieldChangesEne 168 void G4FieldManager::InitialiseFieldChangesEnergy()
189 { 169 {
190 if ( fDetectorField != nullptr ) 170 if ( fDetectorField != nullptr )
191 { 171 {
192 fFieldChangesEnergy = fDetectorField->Doe 172 fFieldChangesEnergy = fDetectorField->DoesFieldChangeEnergy();
193 } 173 }
194 else 174 else
195 { 175 {
196 fFieldChangesEnergy = false; // No fiel 176 fFieldChangesEnergy = false; // No field, no change!
197 } 177 }
198 } 178 }
199 179
200 G4bool G4FieldManager::SetDetectorField(G4Fiel 180 G4bool G4FieldManager::SetDetectorField(G4Field* pDetectorField,
201 G4int 181 G4int failMode )
202 { 182 {
203 G4VIntegrationDriver* driver = nullptr; 183 G4VIntegrationDriver* driver = nullptr;
204 G4EquationOfMotion* equation = nullptr; 184 G4EquationOfMotion* equation = nullptr;
205 // G4bool compatibleField = false; 185 // G4bool compatibleField = false;
206 G4bool ableToSet = false; 186 G4bool ableToSet = false;
207 187
208 fDetectorField = pDetectorField; 188 fDetectorField = pDetectorField;
209 InitialiseFieldChangesEnergy(); 189 InitialiseFieldChangesEnergy();
210 190
211 // Must 'propagate' the field to the depend 191 // Must 'propagate' the field to the dependent classes
212 // 192 //
213 if( fChordFinder != nullptr ) 193 if( fChordFinder != nullptr )
214 { 194 {
215 failMode= std::max( failMode, 1) ; 195 failMode= std::max( failMode, 1) ;
216 // If a chord finder exists, warn in ca 196 // If a chord finder exists, warn in case of error!
217 197
218 driver = fChordFinder->GetIntegrationDriv 198 driver = fChordFinder->GetIntegrationDriver();
219 if( driver != nullptr ) 199 if( driver != nullptr )
220 { 200 {
221 equation = driver->GetEquationOfMotion( 201 equation = driver->GetEquationOfMotion();
222 202
223 // Should check the compatibility betwe 203 // Should check the compatibility between the
224 // field and the equation HERE 204 // field and the equation HERE
225 205
226 if( equation != nullptr ) 206 if( equation != nullptr )
227 { 207 {
228 equation->SetFieldObj(pDetectorField 208 equation->SetFieldObj(pDetectorField);
229 ableToSet = true; 209 ableToSet = true;
230 } 210 }
231 } 211 }
232 } 212 }
233 213
234 if( !ableToSet && (failMode > 0) ) 214 if( !ableToSet && (failMode > 0) )
235 { 215 {
236 // If this fails, report the issue ! 216 // If this fails, report the issue !
237 217
238 G4ExceptionDescription msg; 218 G4ExceptionDescription msg;
239 msg << "Unable to set the field in the de 219 msg << "Unable to set the field in the dependent objects of G4FieldManager"
240 << G4endl; 220 << G4endl;
241 msg << "All the dependent classes must be 221 msg << "All the dependent classes must be fully initialised,"
242 << "before it is possible to call thi 222 << "before it is possible to call this method." << G4endl;
243 msg << "The problem encountered was the f 223 msg << "The problem encountered was the following: " << G4endl;
244 if( fChordFinder == nullptr ) { msg << " 224 if( fChordFinder == nullptr ) { msg << " No ChordFinder. " ; }
245 else if( driver == nullptr ) { msg << " 225 else if( driver == nullptr ) { msg << " No Integration Driver set. ";}
246 else if( equation == nullptr ) { msg << " 226 else if( equation == nullptr ) { msg << " No Equation found. " ; }
247 // else if( !compatibleField ) { msg << " 227 // else if( !compatibleField ) { msg << " Field not compatible. ";}
248 else { msg << " Can NOT find reason for 228 else { msg << " Can NOT find reason for failure. ";}
249 msg << G4endl; 229 msg << G4endl;
250 G4ExceptionSeverity severity = (failMode 230 G4ExceptionSeverity severity = (failMode != 1)
251 ? FatalExcep 231 ? FatalException : JustWarning ;
252 G4Exception("G4FieldManager::SetDetectorF 232 G4Exception("G4FieldManager::SetDetectorField", "Geometry001",
253 severity, msg); 233 severity, msg);
254 } 234 }
255 return ableToSet; 235 return ableToSet;
256 } <<
257 <<
258 G4bool G4FieldManager::SetMaximumEpsilonStep( <<
259 { <<
260 G4bool succeeded= false; <<
261 if( (newEpsMax > 0.0) && ( newEpsMax <= f <<
262 && (fMinAcceptedEpsilon <= newEpsMax ) ) <<
263 { <<
264 if(newEpsMax >= fEpsilonMin) <<
265 { <<
266 fEpsilonMax = newEpsMax; <<
267 succeeded = true; <<
268 if (fVerboseConstruction) <<
269 { <<
270 G4cout << "G4FieldManager/SetEpsMax : <<
271 << " ( Note: unchanged eps_min= <<
272 } <<
273 } <<
274 else <<
275 { <<
276 G4ExceptionDescription erm; <<
277 erm << " Call to set eps_max = " << newE <<
278 << " its value must be at larger or <<
279 erm << " Modifying both to the same valu <<
280 << G4endl <<
281 << " To avoid this warning, please s <<
282 << " 0 < eps_min <= eps_max <= " << <<
283 <<
284 fEpsilonMax = newEpsMax; <<
285 fEpsilonMin = newEpsMax; <<
286 G4String methodName = G4String("G4FieldM <<
287 G4Exception(methodName.c_str(), "Geometr <<
288 } <<
289 } <<
290 else <<
291 { <<
292 G4ExceptionDescription erm; <<
293 G4String paramName("eps_max"); <<
294 ReportBadEpsilonValue(erm, newEpsMax, para <<
295 G4String methodName = G4String("G4FieldMan <<
296 G4Exception(methodName.c_str(), "Geometry0 <<
297 } <<
298 return succeeded; <<
299 } <<
300 <<
301 // ------------------------------------------- <<
302 <<
303 G4bool G4FieldManager::SetMinimumEpsilonStep( <<
304 { <<
305 G4bool succeeded= false; <<
306 <<
307 if( fMinAcceptedEpsilon <= newEpsMin && ne <<
308 { <<
309 fEpsilonMin = newEpsMin; <<
310 //********* <<
311 succeeded= true; <<
312 <<
313 if (fVerboseConstruction) <<
314 { <<
315 G4cout << "G4FieldManager/SetEpsMin : e <<
316 << std::setw(10) << fEpsilonMin < <<
317 } <<
318 if( fEpsilonMax < fEpsilonMin ) <<
319 { <<
320 // Ensure consistency <<
321 G4ExceptionDescription erm; <<
322 erm << "Setting eps_min = " << newEpsMin <<
323 << " For consistency set eps_max= " <<
324 << " ( Old value = " << fEpsilonMax <<
325 fEpsilonMax = fEpsilonMin; <<
326 G4String methodName = G4String("G4FieldM <<
327 G4Exception(methodName.c_str(), "Geometr <<
328 } <<
329 } <<
330 else <<
331 { <<
332 G4ExceptionDescription erm; <<
333 G4String paramName("eps_min"); <<
334 ReportBadEpsilonValue(erm, newEpsMin, para <<
335 G4String methodName = G4String("G4FieldMan <<
336 G4Exception(methodName.c_str(), "Geometry0 <<
337 } <<
338 return succeeded; <<
339 } <<
340 <<
341 // ------------------------------------------- <<
342 <<
343 G4double G4FieldManager::GetMaxAcceptedEpsilon <<
344 { <<
345 return fMaxAcceptedEpsilon; <<
346 } <<
347 <<
348 // ------------------------------------------- <<
349 <<
350 G4bool G4FieldManager::SetMaxAcceptedEpsilon <<
351 // Set value -- within limits <<
352 { <<
353 G4bool success= false; <<
354 // Limit for warning and absolute limit chos <<
355 // investigation of integration with G4Dorma <<
356 if( maxAcceptValue <= fMaxWarningEpsilon ) <<
357 { <<
358 fMaxAcceptedEpsilon= maxAcceptValue; <<
359 success= true; <<
360 } <<
361 else <<
362 { <<
363 G4ExceptionDescription erm; <<
364 G4ExceptionSeverity severity; <<
365 <<
366 G4cout << "G4FieldManager::" << __func__ <<
367 << " Parameters: fMaxAcceptedEpsi <<
368 << " fMaxFinalEpsilon = " << fMaxFi <<
369 <<
370 if( maxAcceptValue <= fMaxFinalEpsilon ) <<
371 { <<
372 success= true; <<
373 fMaxAcceptedEpsilon = maxAcceptValue; <<
374 // Integration is poor, and robustness w <<
375 erm << "Proposed value for maximum-accep <<
376 << " is larger than the recommended <<
377 << G4endl <<
378 << "This may impact the robustness o <<
379 << G4endl <<
380 << "The request was accepted and the <<
381 << " , but future releases are expec <<
382 << " to tighten the limit of accepta <<
383 << fMaxWarningEpsilon << G4endl << G <<
384 << "Suggestion: If you need better p <<
385 << "alternative, low-order RK integr <<
386 << " helix-based methods (for pure B <<
387 << " especially electrons if you nee <<
388 severity= JustWarning; <<
389 } <<
390 else <<
391 { <<
392 fMaxAcceptedEpsilon= fMaxFinalEpsilon; <<
393 erm << " Proposed value for maximum acce <<
394 << " is larger than the top of the r <<
395 << G4endl; <<
396 if( softFailure ) <<
397 { <<
398 erm << " Using the latter value instea <<
399 } <<
400 erm << G4endl; <<
401 erm << " Please adjust to request maxAcc <<
402 << G4endl << G4endl; <<
403 if( !softFailure ) <<
404 { <<
405 erm << " NOTE: you can accept the ceil <<
406 << " warning by using a 2nd argume <<
407 << " in your call to SetMaxAccepte <<
408 } <<
409 severity = softFailure ? JustWarning : F <<
410 // if( softFailure ) severity= JustWarni <<
411 // else severity= FatalException; <<
412 success = false; <<
413 } <<
414 G4String methodName = G4String("G4FieldMan <<
415 G4Exception(methodName.c_str(), "Geometry0 <<
416 } <<
417 return success; <<
418 } <<
419 <<
420 // ------------------------------------------- <<
421 <<
422 void G4FieldManager:: <<
423 ReportBadEpsilonValue(G4ExceptionDescription& <<
424 { <<
425 erm << "Incorrect proposed value of " << nam <<
426 << " Its value is outside the permitted <<
427 << fMinAcceptedEpsilon << " to " << fM <<
428 << " Clarification: " << G4endl; <<
429 G4long oldPrec = erm.precision(); <<
430 if(value < fMinAcceptedEpsilon ) <<
431 { <<
432 erm << " a) The value must be positive an <<
433 << " of the (G4)double type - (" <<
434 << (value < fMinAcceptedEpsilon ? "FA <<
435 << " i.e. std::numeric_limits<G4do <<
436 << std::numeric_limits<G4double>::eps <<
437 << " to ensure that integration " << G <<
438 << " could potentially achieve thi <<
439 << " Minimum accepted eps_min/max <<
440 } <<
441 else if( value > fMaxAcceptedEpsilon) <<
442 { <<
443 erm << " b) It must be smaller than (or e <<
444 << std::setprecision(4) << fMaxAccepte <<
445 << " to ensure robustness of integrati <<
446 << (( value < fMaxAcceptedEpsilon) ? " <<
447 } <<
448 else <<
449 { <<
450 G4bool badRoundoff = (std::fabs(1.0+value) <<
451 erm << " Unknown ERROR case -- extra chec <<
452 erm << " c) as a floating point number (o <<
453 << " ) must be > 1 , (" <<
454 << (badRoundoff ? "FAILED" : "OK" ) < <<
455 << " Now 1+eps_min = " <<
456 << std::setprecision(17) << (1+value) <<
457 << " and (1.0+" << name << ") - <<
458 << std::setprecision(9) << (1.0+value) <<
459 } <<
460 erm.precision(oldPrec); <<
461 } 236 }
462 237