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