| 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 <<
27 // 26 //
28 // Author: John Apostolakis, 10.03.97 - design << 27 // $Id: G4FieldManager.cc,v 1.15 2007/12/07 15:34:10 japost Exp $
>> 28 // GEANT4 tag $Name: geant4-09-03-patch-02 $
>> 29 //
29 // ------------------------------------------- 30 // -------------------------------------------------------------------
30 31
31 #include "G4FieldManager.hh" 32 #include "G4FieldManager.hh"
32 #include "G4Field.hh" 33 #include "G4Field.hh"
33 #include "G4MagneticField.hh" 34 #include "G4MagneticField.hh"
34 #include "G4ChordFinder.hh" 35 #include "G4ChordFinder.hh"
35 #include "G4FieldManagerStore.hh" 36 #include "G4FieldManagerStore.hh"
36 #include "G4SystemOfUnits.hh" <<
37 37
38 G4ThreadLocal G4FieldManager* G4FieldManager:: << 38 G4FieldManager::G4FieldManager(G4Field *detectorField,
39 G4double G4FieldManager::fDefault_Delta_One_St << 39 G4ChordFinder *pChordFinder,
40 G4double G4FieldManager::fDefault_Delta_Inters << 40 G4bool fieldChangesEnergy
41 G4bool G4FieldManager::fVerboseConstruction= << 41 )
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 <<
50 G4FieldManager::G4FieldManager(G4Field* detect <<
51 G4ChordFinder* <<
52 G4bool fieldCha <<
53 : fDetectorField(detectorField), 42 : fDetectorField(detectorField),
54 fChordFinder(pChordFinder), 43 fChordFinder(pChordFinder),
55 fDelta_One_Step_Value( fDefault_Delta_One << 44 fAllocatedChordFinder(false),
56 fDelta_Intersection_Val( fDefault_Delta_I << 45 fDefault_Delta_One_Step_Value(0.01*mm),
>> 46 fDefault_Delta_Intersection_Val(0.001*mm),
>> 47 fEpsilonMinDefault(5.0e-5),
>> 48 fEpsilonMaxDefault(0.001),
57 fEpsilonMin( fEpsilonMinDefault ), 49 fEpsilonMin( fEpsilonMinDefault ),
58 fEpsilonMax( fEpsilonMaxDefault) 50 fEpsilonMax( fEpsilonMaxDefault)
59 { 51 {
60 if ( detectorField != nullptr ) << 52 fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
61 { << 53 fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
62 fFieldChangesEnergy = detectorField->Does << 54 if ( detectorField )
63 } << 55 fFieldChangesEnergy= detectorField->DoesFieldChangeEnergy();
64 else 56 else
65 { << 57 fFieldChangesEnergy= fieldChangesEnergy;
66 fFieldChangesEnergy = fieldChangesEnergy; <<
67 } <<
68 <<
69 if( fVerboseConstruction) <<
70 { <<
71 G4cout << "G4FieldManager/ctor#1 fEpsilon <<
72 } <<
73 58
74 // Add to store 59 // Add to store
75 // <<
76 G4FieldManagerStore::Register(this); 60 G4FieldManagerStore::Register(this);
>> 61
77 } 62 }
78 63
79 G4FieldManager::G4FieldManager(G4MagneticField << 64 G4FieldManager::G4FieldManager(G4MagneticField *detectorField)
80 : fDetectorField(detectorField), fAllocated 65 : fDetectorField(detectorField), fAllocatedChordFinder(true),
81 fDelta_One_Step_Value( fDefault_Delta_O << 66 fFieldChangesEnergy(false),
82 fDelta_Intersection_Val( fDefault_Delta_I << 67 fDefault_Delta_One_Step_Value(0.01*mm),
>> 68 fDefault_Delta_Intersection_Val(0.001*mm),
>> 69 fEpsilonMinDefault(5.0e-5),
>> 70 fEpsilonMaxDefault(0.001),
83 fEpsilonMin( fEpsilonMinDefault ), 71 fEpsilonMin( fEpsilonMinDefault ),
84 fEpsilonMax( fEpsilonMaxDefault ) << 72 fEpsilonMax( fEpsilonMaxDefault)
85 { 73 {
86 fChordFinder = new G4ChordFinder( detectorF << 74 fChordFinder= new G4ChordFinder( detectorField );
87 << 75 fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
88 if( fVerboseConstruction ) << 76 fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
89 { <<
90 G4cout << "G4FieldManager/ctor#2 fEpsilon <<
91 } <<
92 // Add to store 77 // Add to store
93 // <<
94 G4FieldManagerStore::Register(this); 78 G4FieldManagerStore::Register(this);
95 } 79 }
96 80
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 th <<
109 // any chordfinder now <<
110 // <<
111 aFM = new G4FieldManager( aField , nul <<
112 <<
113 // Check if originally we have the fAl <<
114 // set, in case, call chord constructo <<
115 // <<
116 if ( fAllocatedChordFinder ) <<
117 { <<
118 aFM->CreateChordFinder( dynamic_ca <<
119 } <<
120 else <<
121 { <<
122 // Chord was specified by user, sh <<
123 // TODO: For the moment copy point <<
124 // if cloning of ChordFinder <<
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_ <<
135 aFM->fDelta_One_Step_Value = fDelta_On <<
136 // TODO: Should we really add to the <<
137 // Who will use this? <<
138 } <<
139 catch ( ... ) <<
140 { <<
141 // Failed creating clone: probably use <<
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 M <<
151 return aFM; <<
152 } <<
153 <<
154 void G4FieldManager::ConfigureForTrack( const 81 void G4FieldManager::ConfigureForTrack( const G4Track * )
155 { 82 {
156 // Default is to do nothing! 83 // Default is to do nothing!
>> 84 ;
157 } 85 }
158 86
159 G4FieldManager::~G4FieldManager() 87 G4FieldManager::~G4FieldManager()
160 { 88 {
161 if( fAllocatedChordFinder ) << 89 if( fAllocatedChordFinder ){
162 { <<
163 delete fChordFinder; 90 delete fChordFinder;
164 } 91 }
165 G4FieldManagerStore::DeRegister(this); 92 G4FieldManagerStore::DeRegister(this);
166 } 93 }
167 94
168 void 95 void
169 G4FieldManager::CreateChordFinder(G4MagneticFi << 96 G4FieldManager::CreateChordFinder(G4MagneticField *detectorMagField)
170 { 97 {
171 if ( fAllocatedChordFinder ) 98 if ( fAllocatedChordFinder )
172 { <<
173 delete fChordFinder; 99 delete fChordFinder;
174 } << 100 fChordFinder= new G4ChordFinder( detectorMagField );
175 fAllocatedChordFinder = false; << 101 fAllocatedChordFinder= true;
176 <<
177 if( detectorMagField != nullptr ) <<
178 { <<
179 fChordFinder = new G4ChordFinder( detect <<
180 fAllocatedChordFinder = true; <<
181 } <<
182 else <<
183 { <<
184 fChordFinder = nullptr; <<
185 } <<
186 } <<
187 <<
188 void G4FieldManager::InitialiseFieldChangesEne <<
189 { <<
190 if ( fDetectorField != nullptr ) <<
191 { <<
192 fFieldChangesEnergy = fDetectorField->Doe <<
193 } <<
194 else <<
195 { <<
196 fFieldChangesEnergy = false; // No fiel <<
197 } <<
198 } <<
199 <<
200 G4bool G4FieldManager::SetDetectorField(G4Fiel <<
201 G4int <<
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 depend <<
212 // <<
213 if( fChordFinder != nullptr ) <<
214 { <<
215 failMode= std::max( failMode, 1) ; <<
216 // If a chord finder exists, warn in ca <<
217 <<
218 driver = fChordFinder->GetIntegrationDriv <<
219 if( driver != nullptr ) <<
220 { <<
221 equation = driver->GetEquationOfMotion( <<
222 <<
223 // Should check the compatibility betwe <<
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 de <<
240 << G4endl; <<
241 msg << "All the dependent classes must be <<
242 << "before it is possible to call thi <<
243 msg << "The problem encountered was the f <<
244 if( fChordFinder == nullptr ) { msg << " <<
245 else if( driver == nullptr ) { msg << " <<
246 else if( equation == nullptr ) { msg << " <<
247 // else if( !compatibleField ) { msg << " <<
248 else { msg << " Can NOT find reason for <<
249 msg << G4endl; <<
250 G4ExceptionSeverity severity = (failMode <<
251 ? FatalExcep <<
252 G4Exception("G4FieldManager::SetDetectorF <<
253 severity, msg); <<
254 } <<
255 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 } 102 }
340 103
341 // ------------------------------------------- << 104 G4bool G4FieldManager::SetDetectorField(G4Field *pDetectorField)
342 <<
343 G4double G4FieldManager::GetMaxAcceptedEpsilon <<
344 { 105 {
345 return fMaxAcceptedEpsilon; << 106 fDetectorField= pDetectorField;
346 } <<
347 107
348 // ------------------------------------------- << 108 if ( pDetectorField )
>> 109 fFieldChangesEnergy= pDetectorField->DoesFieldChangeEnergy();
>> 110 else
>> 111 fFieldChangesEnergy= false; // No field
349 112
350 G4bool G4FieldManager::SetMaxAcceptedEpsilon << 113 return false;
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 } 114 }
419 115
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 } <<
462 116