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-01 $ >> 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