Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // 32 // File name: G4VMultipleScattering 33 // 34 // Author: Vladimir Ivanchenko on base 35 // 36 // Creation date: 25.03.2003 37 // 38 // Modifications: 39 // 40 // 16-07-03 Use G4VMscModel interface (V.Ivanc 41 // 03-11-03 Fix initialisation problem in Retr 42 // 04-11-03 Update PrintInfoDefinition (V.Ivan 43 // 01-03-04 SampleCosineTheta signature change 44 // 22-04-04 SampleCosineTheta signature change 45 // 27-08-04 Add InitialiseForRun method (V.Iva 46 // 08-11-04 Migration to new interface of Stor 47 // 11-03-05 Shift verbose level by 1 (V.Ivantc 48 // 15-04-05 optimize internal interface (V.Iva 49 // 15-04-05 remove boundary flag (V.Ivanchenko 50 // 27-10-05 introduce virtual function MscStep 51 // 12-04-07 Add verbosity at destruction (V.Iv 52 // 27-10-07 Virtual functions moved to source 53 // 11-03-08 Set skin value does not effect ste 54 // 24-06-09 Removed hidden bin in G4PhysicsVec 55 // 04-06-13 Adoptation to MT mode (V.Ivanchenk 56 // 57 58 // ------------------------------------------- 59 // 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oo 62 63 #include "G4VMultipleScattering.hh" 64 #include "G4PhysicalConstants.hh" 65 #include "G4SystemOfUnits.hh" 66 #include "G4LossTableManager.hh" 67 #include "G4MaterialCutsCouple.hh" 68 #include "G4Step.hh" 69 #include "G4ParticleDefinition.hh" 70 #include "G4VEmFluctuationModel.hh" 71 #include "G4UnitsTable.hh" 72 #include "G4ProductionCutsTable.hh" 73 #include "G4Electron.hh" 74 #include "G4GenericIon.hh" 75 #include "G4TransportationManager.hh" 76 #include "G4SafetyHelper.hh" 77 #include "G4ParticleTable.hh" 78 #include "G4ProcessVector.hh" 79 #include "G4ProcessManager.hh" 80 #include "G4LossTableBuilder.hh" 81 #include "G4EmTableUtil.hh" 82 #include <iostream> 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 85 86 G4VMultipleScattering::G4VMultipleScattering(c 87 : G4VContinuousDiscreteProcess("msc", fElect 88 fNewPosition(0.,0.,0.), 89 fNewDirection(0.,0.,1.) 90 { 91 theParameters = G4EmParameters::Instance(); 92 SetVerboseLevel(1); 93 SetProcessSubType(fMultipleScattering); 94 95 lowestKinEnergy = 10*CLHEP::eV; 96 97 geomMin = 0.05*CLHEP::nm; 98 minDisplacement2 = geomMin*geomMin; 99 100 pParticleChange = &fParticleChange; 101 102 modelManager = new G4EmModelManager(); 103 emManager = G4LossTableManager::Instance(); 104 mscModels.reserve(2); 105 emManager->Register(this); 106 } 107 108 //....oooOO0OOooo........oooOO0OOooo........oo 109 110 G4VMultipleScattering::~G4VMultipleScattering( 111 { 112 delete modelManager; 113 emManager->DeRegister(this); 114 } 115 116 //....oooOO0OOooo........oooOO0OOooo........oo 117 118 void G4VMultipleScattering::AddEmModel(G4int o 119 const G 120 { 121 if(nullptr == ptr) { return; } 122 G4VEmFluctuationModel* fm = nullptr; 123 modelManager->AddEmModel(order, ptr, fm, reg 124 ptr->SetParticleChange(pParticleChange); 125 } 126 127 //....oooOO0OOooo........oooOO0OOooo........oo 128 129 void G4VMultipleScattering::SetEmModel(G4VMscM 130 { 131 if(nullptr == ptr) { return; } 132 if(!mscModels.empty()) { 133 for(auto & msc : mscModels) { if(msc == pt 134 } 135 mscModels.push_back(ptr); 136 } 137 138 //....oooOO0OOooo........oooOO0OOooo........oo 139 140 void 141 G4VMultipleScattering::PreparePhysicsTable(con 142 { 143 G4bool master = emManager->IsMaster(); 144 if (nullptr == firstParticle) { firstParticl 145 146 emManager->PreparePhysicsTable(&part, this); 147 currParticle = nullptr; 148 149 if(firstParticle == &part) { 150 baseMat = emManager->GetTableBuilder()->Ge 151 G4EmTableUtil::PrepareMscProcess(this, par 152 stepLimit, facrange, 153 latDisplacement, master, 154 isIon, baseMat); 155 156 numberOfModels = modelManager->NumberOfMod 157 currentModel = GetModelByIndex(0); 158 159 if (nullptr == safetyHelper) { 160 safetyHelper = G4TransportationManager:: 161 ->GetSafetyHelper(); 162 safetyHelper->InitialiseHelper(); 163 } 164 } 165 } 166 167 //....oooOO0OOooo........oooOO0OOooo........oo 168 169 void G4VMultipleScattering::BuildPhysicsTable( 170 { 171 G4bool master = emManager->IsMaster(); 172 173 if(firstParticle == &part) { 174 emManager->BuildPhysicsTable(&part); 175 } 176 const G4VMultipleScattering* ptr = this; 177 if(!master) { 178 ptr = static_cast<const G4VMultipleScatter 179 } 180 181 G4EmTableUtil::BuildMscProcess(this, ptr, pa 182 numberOfModels, master); 183 } 184 185 //....oooOO0OOooo........oooOO0OOooo........oo 186 187 void G4VMultipleScattering::StreamInfo(std::os 188 const G4ParticleDefinition& 189 { 190 G4String indent = (rst ? " " : ""); 191 outFile << G4endl << indent << GetProcessNam 192 if (!rst) outFile << " for " << part.GetPart 193 outFile << " SubType= " << GetProcessSubTy 194 modelManager->DumpModelList(outFile, verbose 195 } 196 197 //....oooOO0OOooo........oooOO0OOooo........oo 198 199 void G4VMultipleScattering::StartTracking(G4Tr 200 { 201 G4VEnergyLossProcess* eloss = nullptr; 202 if(track->GetParticleDefinition() != currPar 203 currParticle = track->GetParticleDefinitio 204 fIonisation = emManager->GetEnergyLossProc 205 eloss = fIonisation; 206 } 207 for(G4int i=0; i<numberOfModels; ++i) { 208 G4VMscModel* msc = GetModelByIndex(i); 209 msc->StartTracking(track); 210 if(nullptr != eloss) { 211 msc->SetIonisation(eloss, currParticle); 212 } 213 } 214 } 215 216 //....oooOO0OOooo........oooOO0OOooo........oo 217 218 G4double G4VMultipleScattering::AlongStepGetPh 219 const G4Track& tr 220 G4double, 221 G4double currentM 222 G4double&, 223 G4GPILSelection* 224 { 225 // get Step limit proposed by the process 226 *selection = NotCandidateForSelection; 227 physStepLimit = gPathLength = tPathLength = 228 229 G4double ekin = track.GetKineticEnergy(); 230 /* 231 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 232 << " " << currParticle->GetParticleN 233 << " currMod " << currentModel 234 << G4endl; 235 */ 236 // isIon flag is used only to select a model 237 if(isIon) { 238 ekin *= proton_mass_c2/track.GetParticleDe 239 } 240 const G4MaterialCutsCouple* couple = track.G 241 242 // select new model, static cast is possible 243 if(1 < numberOfModels) { 244 currentModel = 245 static_cast<G4VMscModel*>(SelectModel(ek 246 } 247 currentModel->SetCurrentCouple(couple); 248 // msc is active is model is active, energy 249 // and step size is above the limit; 250 // if it is active msc may limit the step 251 if(currentModel->IsActive(ekin) && tPathLeng 252 && ekin >= lowestKinEnergy) { 253 isActive = true; 254 tPathLength = 255 currentModel->ComputeTruePathLengthLimit 256 if (tPathLength < physStepLimit) { 257 *selection = CandidateForSelection; 258 } 259 } else { 260 isActive = false; 261 gPathLength = DBL_MAX; 262 } 263 264 //if(currParticle->GetPDGMass() > GeV) 265 /* 266 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 267 << " " << currParticle->GetParticleN 268 << " gPathLength= " << gPathLength 269 << " tPathLength= " << tPathLength 270 << " currentMinimalStep= " << current 271 << " isActive " << isActive << G4endl 272 */ 273 return gPathLength; 274 } 275 276 //....oooOO0OOooo........oooOO0OOooo........oo 277 278 G4double 279 G4VMultipleScattering::PostStepGetPhysicalInte 280 const G4Track&, G4double, G4Forc 281 { 282 *condition = NotForced; 283 return DBL_MAX; 284 } 285 286 //....oooOO0OOooo........oooOO0OOooo........oo 287 288 G4VParticleChange* 289 G4VMultipleScattering::AlongStepDoIt(const G4T 290 { 291 fParticleChange.InitialiseMSC(track, step); 292 fNewPosition = fParticleChange.GetProposedPo 293 fPositionChanged = false; 294 295 G4double geomLength = step.GetStepLength(); 296 297 // very small step - no msc 298 if(!isActive) { 299 tPathLength = geomLength; 300 301 // sample msc 302 } else { 303 G4double range = 304 currentModel->GetRange(currParticle,trac 305 track.GetMaterial 306 307 tPathLength = currentModel->ComputeTrueSte 308 309 /* 310 if(currParticle->GetPDGMass() > 0.9*GeV) 311 G4cout << "G4VMsc::AlongStepDoIt: GeomLeng 312 << geomLength 313 << " tPathLength= " << tPathLength 314 << " physStepLimit= " << physStepLi 315 << " dr= " << range - tPathLength 316 << " ekin= " << track.GetKineticEne 317 */ 318 // protection against wrong t->g->t conver 319 tPathLength = std::min(tPathLength, physSt 320 321 // do not sample scattering at the last or 322 if(tPathLength < range && tPathLength > ge 323 324 static const G4double minSafety = 1.20*C 325 static const G4double sFact = 0.99; 326 327 G4ThreeVector displacement = currentMode 328 step.GetPostStepPoint()->GetMomentumDi 329 330 G4double r2 = displacement.mag2(); 331 //G4cout << " R= " << sqrt(r2) << " R 332 // << " flag= " << fDispBeyondSafety 333 if(r2 > minDisplacement2) { 334 335 fPositionChanged = true; 336 G4double dispR = std::sqrt(r2); 337 G4double postSafety = 338 sFact*safetyHelper->ComputeSafety(fN 339 //G4cout<<" R= "<< dispR<<" postSaf 340 341 // far away from geometry boundary 342 if(postSafety > 0.0 && dispR <= postSa 343 fNewPosition += displacement; 344 345 //near the boundary 346 } else { 347 // displaced point is definitely wit 348 //G4cout<<" R= "<<dispR<<" postSa 349 if(dispR < postSafety) { 350 fNewPosition += displacement; 351 352 // reduced displacement 353 } else if(postSafety > geomMin) { 354 fNewPosition += displacement*(post 355 356 // very small postSafety 357 } else { 358 fPositionChanged = false; 359 } 360 } 361 if(fPositionChanged) { 362 safetyHelper->ReLocateWithinVolume(f 363 fParticleChange.ProposePosition(fNew 364 } 365 } 366 } 367 } 368 fParticleChange.ProposeTrueStepLength(tPathL 369 return &fParticleChange; 370 } 371 372 //....oooOO0OOooo........oooOO0OOooo........oo 373 374 G4double G4VMultipleScattering::GetContinuousS 375 const G 376 G4doubl 377 G4doubl 378 G4doubl 379 { 380 G4GPILSelection selection = NotCandidateForS 381 G4double x = AlongStepGetPhysicalInteraction 382 383 384 385 return x; 386 } 387 388 //....oooOO0OOooo........oooOO0OOooo........oo 389 390 G4double G4VMultipleScattering::ContinuousStep 391 const G 392 G4doubl 393 G4doubl 394 G4doubl 395 { 396 return GetContinuousStepLimit(track,previous 397 currentSafety) 398 } 399 400 //....oooOO0OOooo........oooOO0OOooo........oo 401 402 G4double G4VMultipleScattering::GetMeanFreePat 403 const G4Track&, G4double, G4Forc 404 { 405 *condition = Forced; 406 return DBL_MAX; 407 } 408 409 //....oooOO0OOooo........oooOO0OOooo........oo 410 411 G4bool 412 G4VMultipleScattering::StorePhysicsTable(const 413 const 414 G4boo 415 { 416 G4bool yes = true; 417 if(part != firstParticle || !emManager->IsMa 418 419 return G4EmTableUtil::StoreMscTable(this, pa 420 numberOfModels, verboseLevel, 421 ascii); 422 } 423 424 //....oooOO0OOooo........oooOO0OOooo........oo 425 426 G4bool 427 G4VMultipleScattering::RetrievePhysicsTable(co 428 co 429 G4 430 { 431 return true; 432 } 433 434 //....oooOO0OOooo........oooOO0OOooo........oo 435 436 void G4VMultipleScattering::ProcessDescription 437 { 438 if(nullptr != firstParticle) { 439 StreamInfo(outFile, *firstParticle, true); 440 } 441 } 442 443 //....oooOO0OOooo........oooOO0OOooo........oo 444 445