| 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 107059 2017-11-01 14:58:16Z gcosmo $
>> 28 //
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" <<
37 36
38 G4ThreadLocal G4FieldManager* G4FieldManager:: << 37 G4FieldManager::G4FieldManager(G4Field *detectorField,
39 G4double G4FieldManager::fDefault_Delta_One_St << 38 G4ChordFinder *pChordFinder,
40 G4double G4FieldManager::fDefault_Delta_Inters << 39 G4bool fieldChangesEnergy
41 G4bool G4FieldManager::fVerboseConstruction= << 40 )
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), 41 : fDetectorField(detectorField),
54 fChordFinder(pChordFinder), 42 fChordFinder(pChordFinder),
55 fDelta_One_Step_Value( fDefault_Delta_One << 43 fAllocatedChordFinder(false),
56 fDelta_Intersection_Val( fDefault_Delta_I << 44 fEpsilonMinDefault(5.0e-5),
>> 45 fEpsilonMaxDefault(0.001),
>> 46 fDefault_Delta_One_Step_Value(0.01), // mm
>> 47 fDefault_Delta_Intersection_Val(0.001), // mm
57 fEpsilonMin( fEpsilonMinDefault ), 48 fEpsilonMin( fEpsilonMinDefault ),
58 fEpsilonMax( fEpsilonMaxDefault) 49 fEpsilonMax( fEpsilonMaxDefault)
59 { 50 {
60 if ( detectorField != nullptr ) << 51 fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
61 { << 52 fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
62 fFieldChangesEnergy = detectorField->Does << 53 if ( detectorField )
63 } << 54 fFieldChangesEnergy= detectorField->DoesFieldChangeEnergy();
64 else 55 else
65 { << 56 fFieldChangesEnergy= fieldChangesEnergy;
66 fFieldChangesEnergy = fieldChangesEnergy; <<
67 } <<
68 <<
69 if( fVerboseConstruction) <<
70 { <<
71 G4cout << "G4FieldManager/ctor#1 fEpsilon <<
72 } <<
73 57
74 // Add to store 58 // Add to store
75 // <<
76 G4FieldManagerStore::Register(this); 59 G4FieldManagerStore::Register(this);
77 } 60 }
78 61
79 G4FieldManager::G4FieldManager(G4MagneticField << 62 G4FieldManager::G4FieldManager(G4MagneticField *detectorField)
80 : fDetectorField(detectorField), fAllocated 63 : fDetectorField(detectorField), fAllocatedChordFinder(true),
81 fDelta_One_Step_Value( fDefault_Delta_O << 64 fEpsilonMinDefault(5.0e-5),
82 fDelta_Intersection_Val( fDefault_Delta_I << 65 fEpsilonMaxDefault(0.001),
>> 66 fFieldChangesEnergy(false),
>> 67 fDefault_Delta_One_Step_Value(0.01), // mm
>> 68 fDefault_Delta_Intersection_Val(0.001), // mm
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 fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
88 if( fVerboseConstruction ) << 74 fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
89 { <<
90 G4cout << "G4FieldManager/ctor#2 fEpsilon <<
91 } <<
92 // Add to store 75 // Add to store
93 // <<
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 = 0;
100 G4FieldManager* aFM = nullptr; << 82 G4FieldManager* aFM = 0;
101 G4ChordFinder* aCF = nullptr; << 83 G4ChordFinder* aCF = 0;
102 try { 84 try {
103 if ( fDetectorField != nullptr ) << 85 if ( this->fDetectorField )
104 { << 86 aField = this->fDetectorField->Clone();
105 aField = fDetectorField->Clone(); << 87
106 } << 88 //Create a new field manager, note that we do not set any chordfinder now.
>> 89 aFM = new G4FieldManager( aField , 0 , this->fFieldChangesEnergy );
107 90
108 // Create a new field manager, note th << 91 //Check if orignally we have the fAllocatedChordFinder variable set, in case, call chord
109 // any chordfinder now << 92 //constructor
110 // << 93 if ( this->fAllocatedChordFinder )
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 { 94 {
118 aFM->CreateChordFinder( dynamic_ca 95 aFM->CreateChordFinder( dynamic_cast<G4MagneticField*>(aField) );
119 } 96 }
120 else 97 else
121 { 98 {
122 // Chord was specified by user, sh << 99 //Chord was specified by user, should we clone?
123 // TODO: For the moment copy point << 100 //TODO: For the moment copy pointer, to be understood if cloning of ChordFinder is needed
124 // if cloning of ChordFinder << 101 aCF = this->fChordFinder;/*->Clone*/
125 // <<
126 aCF = fChordFinder; /*->Clone*/ <<
127 aFM->fChordFinder = aCF; 102 aFM->fChordFinder = aCF;
128 } 103 }
129 << 104 //Copy values of other variables
130 // Copy values of other variables << 105 aFM->fEpsilonMax = this->fEpsilonMax;
131 << 106 aFM->fEpsilonMin = this->fEpsilonMin;
132 aFM->fEpsilonMax = fEpsilonMax; << 107 aFM->fDefault_Delta_Intersection_Val = this->fDefault_Delta_Intersection_Val;
133 aFM->fEpsilonMin = fEpsilonMin; << 108 aFM->fDefault_Delta_One_Step_Value = this->fDefault_Delta_One_Step_Value;
134 aFM->fDelta_Intersection_Val = fDelta_ << 109 aFM->fDelta_Intersection_Val = this->fDelta_Intersection_Val;
135 aFM->fDelta_One_Step_Value = fDelta_On << 110 aFM->fDelta_One_Step_Value = this->fDelta_One_Step_Value;
136 // TODO: Should we really add to the << 111 //TODO: Should we really add to the store the cloned FM? Who will use this?
137 // Who will use this? <<
138 } 112 }
139 catch ( ... ) 113 catch ( ... )
140 { 114 {
141 // Failed creating clone: probably use << 115 //Failed creating clone: probably user did not implement Clone method
142 // in derived classes? << 116 //in derived classes?
143 // Perform clean-up after ourselves... << 117 //Perform clean-up after ourselves...
144 delete aField; 118 delete aField;
145 delete aFM; 119 delete aFM;
146 delete aCF; 120 delete aCF;
147 throw; 121 throw;
148 } 122 }
149 <<
150 G4cout << "G4FieldManager/clone fEpsilon M <<
151 return aFM; 123 return aFM;
152 } 124 }
153 125
154 void G4FieldManager::ConfigureForTrack( const 126 void G4FieldManager::ConfigureForTrack( const G4Track * )
155 { 127 {
156 // Default is to do nothing! 128 // Default is to do nothing!
>> 129 ;
157 } 130 }
158 131
159 G4FieldManager::~G4FieldManager() 132 G4FieldManager::~G4FieldManager()
160 { 133 {
161 if( fAllocatedChordFinder ) << 134 if( fAllocatedChordFinder ){
162 { <<
163 delete fChordFinder; 135 delete fChordFinder;
164 } 136 }
165 G4FieldManagerStore::DeRegister(this); 137 G4FieldManagerStore::DeRegister(this);
166 } 138 }
167 139
168 void 140 void
169 G4FieldManager::CreateChordFinder(G4MagneticFi << 141 G4FieldManager::CreateChordFinder(G4MagneticField *detectorMagField)
170 { 142 {
171 if ( fAllocatedChordFinder ) 143 if ( fAllocatedChordFinder )
172 { <<
173 delete fChordFinder; 144 delete fChordFinder;
174 } << 145 fChordFinder= new G4ChordFinder( detectorMagField );
175 fAllocatedChordFinder = false; << 146 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 } 147 }
187 148
188 void G4FieldManager::InitialiseFieldChangesEne 149 void G4FieldManager::InitialiseFieldChangesEnergy()
189 { 150 {
190 if ( fDetectorField != nullptr ) << 151 if ( fDetectorField )
191 { << 152 fFieldChangesEnergy= fDetectorField->DoesFieldChangeEnergy();
192 fFieldChangesEnergy = fDetectorField->Doe <<
193 } <<
194 else 153 else
195 { << 154 fFieldChangesEnergy= false; // No field , no change!
196 fFieldChangesEnergy = false; // No fiel <<
197 } <<
198 } 155 }
199 156
200 G4bool G4FieldManager::SetDetectorField(G4Fiel << 157 G4bool G4FieldManager::SetDetectorField(G4Field *pDetectorField, int failMode )
201 G4int <<
202 { 158 {
203 G4VIntegrationDriver* driver = nullptr; 159 G4VIntegrationDriver* driver = nullptr;
204 G4EquationOfMotion* equation = nullptr; 160 G4EquationOfMotion* equation = nullptr;
205 // G4bool compatibleField = false; << 161 // G4bool compatibleField= false;
206 G4bool ableToSet = false; << 162 G4bool ableToSet= false;
207 163
208 fDetectorField = pDetectorField; << 164 fDetectorField= pDetectorField;
209 InitialiseFieldChangesEnergy(); 165 InitialiseFieldChangesEnergy();
210 166
211 // Must 'propagate' the field to the depend 167 // Must 'propagate' the field to the dependent classes
212 // 168 //
213 if( fChordFinder != nullptr ) << 169 if( fChordFinder )
214 { 170 {
215 failMode= std::max( failMode, 1) ; << 171 failMode= std::max( failMode, 1) ; // If a chord finder exists, warn in case of error!
216 // If a chord finder exists, warn in ca <<
217 172
218 driver = fChordFinder->GetIntegrationDriv 173 driver = fChordFinder->GetIntegrationDriver();
219 if( driver != nullptr ) << 174 if( driver ){
220 { << 175 equation = driver->GetEquationOfMotion();
221 equation = driver->GetEquationOfMotion( << 176 // Should check the compatibility between the field and the equation HERE
222 << 177 if( equation ) {
223 // Should check the compatibility betwe << 178 equation->SetFieldObj(pDetectorField);
224 // field and the equation HERE << 179 ableToSet = true;
225 << 180 }
226 if( equation != nullptr ) <<
227 { <<
228 equation->SetFieldObj(pDetectorField <<
229 ableToSet = true; <<
230 } <<
231 } 181 }
232 } 182 }
233 183
234 if( !ableToSet && (failMode > 0) ) 184 if( !ableToSet && (failMode > 0) )
235 { 185 {
236 // If this fails, report the issue ! << 186 // If this fails, report the issue !
237 << 187 G4ExceptionDescription msg;
238 G4ExceptionDescription msg; << 188 msg << "Unable to set the field in the dependent objects of G4FieldManager" << G4endl;
239 msg << "Unable to set the field in the de << 189 msg << "All the dependent classes must be fully initialised, before it is possible to call this method." << G4endl;
240 << G4endl; << 190 msg << "The problem encountered was the following: " << G4endl;
241 msg << "All the dependent classes must be << 191 if( !fChordFinder ) { msg << " No ChordFinder. " ; }
242 << "before it is possible to call thi << 192 else if( !driver) { msg << " No Integration Driver set. ";}
243 msg << "The problem encountered was the f << 193 else if( !equation) { msg << " No Equation found. " ; }
244 if( fChordFinder == nullptr ) { msg << " << 194 // else if( !compatibleField ) { msg << " Field not compatible. ";}
245 else if( driver == nullptr ) { msg << " << 195 else { msg << " Can NOT find reason for failure. ";}
246 else if( equation == nullptr ) { msg << " << 196 msg << G4endl;
247 // else if( !compatibleField ) { msg << " << 197 G4ExceptionSeverity severity= (failMode != 1) ? FatalException : JustWarning ;
248 else { msg << " Can NOT find reason for << 198 G4Exception("G4FieldManager::SetDetectorField", "Geometry001",
249 msg << G4endl; << 199 severity, msg);
250 G4ExceptionSeverity severity = (failMode <<
251 ? FatalExcep <<
252 G4Exception("G4FieldManager::SetDetectorF <<
253 severity, msg); <<
254 } 200 }
255 return ableToSet; 201 return ableToSet;
256 } 202 }
257 203
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 } <<
462 204