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 // G4TransportationWithMsc 27 // 28 // Class Description: 29 // 30 // It is a generic process of transportation with multiple scattering included 31 // in the step limitation and propagation. 32 // 33 // Original author: Jonas Hahnfeld, 2022 34 35 // ------------------------------------------------------------------- 36 // 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 40 #include "G4TransportationWithMsc.hh" 41 42 #include "G4DynamicParticle.hh" 43 #include "G4Electron.hh" 44 #include "G4EmConfigurator.hh" 45 #include "G4EmDataHandler.hh" 46 #include "G4LossTableBuilder.hh" 47 #include "G4LossTableManager.hh" 48 #include "G4ParticleChangeForGamma.hh" 49 #include "G4ParticleChangeForMSC.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicsTableHelper.hh" 52 #include "G4PhysicsVector.hh" 53 #include "G4ProductionCutsTable.hh" 54 #include "G4Step.hh" 55 #include "G4StepPoint.hh" 56 #include "G4StepStatus.hh" 57 #include "G4Track.hh" 58 #include "G4TransportationProcessType.hh" 59 #include "G4VMscModel.hh" 60 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 63 static constexpr G4double kLowestKinEnergy = 10 * CLHEP::eV; 64 static constexpr G4double kGeomMin = 0.05 * CLHEP::nm; 65 static constexpr G4double kMinDisplacement2 = kGeomMin * kGeomMin; 66 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 68 69 G4TransportationWithMsc::G4TransportationWithMsc(ScatteringType type, G4int verbosity) 70 : G4Transportation(verbosity, "TransportationWithMsc"), fType(type) 71 { 72 SetProcessSubType(static_cast<G4int>(TRANSPORTATION_WITH_MSC)); 73 SetVerboseLevel(1); 74 75 fEmManager = G4LossTableManager::Instance(); 76 fModelManager = new G4EmModelManager; 77 78 if (type == ScatteringType::MultipleScattering) { 79 fParticleChangeForMSC = new G4ParticleChangeForMSC; 80 } 81 else if (type == ScatteringType::SingleScattering) { 82 fParticleChangeForSS = new G4ParticleChangeForGamma; 83 fSecondariesSS = new std::vector<G4DynamicParticle*>; 84 } 85 86 G4ThreeVector zero; 87 fSubStepDynamicParticle = new G4DynamicParticle(G4Electron::Definition(), zero); 88 fSubStepTrack = new G4Track(fSubStepDynamicParticle, 0, zero); 89 fSubStep = new G4Step; 90 fSubStepTrack->SetStep(fSubStep); 91 } 92 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 94 95 G4TransportationWithMsc::~G4TransportationWithMsc() 96 { 97 delete fModelManager; 98 delete fParticleChangeForMSC; 99 delete fEmData; 100 delete fParticleChangeForSS; 101 delete fSecondariesSS; 102 103 // fSubStepDynamicParticle is owned and also deleted by fSubStepTrack! 104 delete fSubStepTrack; 105 delete fSubStep; 106 } 107 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 109 110 void G4TransportationWithMsc::AddMscModel(G4VMscModel* mscModel, G4int order, 111 const G4Region* region) 112 { 113 if (fType != ScatteringType::MultipleScattering) { 114 G4Exception("G4TransportationWithMsc::AddMscModel", "em0051", FatalException, 115 "not allowed unless type == MultipleScattering"); 116 } 117 118 fModelManager->AddEmModel(order, mscModel, nullptr, region); 119 mscModel->SetParticleChange(fParticleChangeForMSC); 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 124 void G4TransportationWithMsc::AddSSModel(G4VEmModel* model, G4int order, const G4Region* region) 125 { 126 if (fType != ScatteringType::SingleScattering) { 127 G4Exception("G4TransportationWithMsc::AddSSModel", "em0051", FatalException, 128 "not allowed unless type == SingleScattering"); 129 } 130 131 fModelManager->AddEmModel(order, model, nullptr, region); 132 model->SetPolarAngleLimit(0.0); 133 model->SetParticleChange(fParticleChangeForSS); 134 } 135 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 137 138 void G4TransportationWithMsc::PreparePhysicsTable(const G4ParticleDefinition& part) 139 { 140 if (nullptr == fFirstParticle) { 141 fFirstParticle = ∂ 142 G4VMultipleScattering* ptr = nullptr; 143 auto emConfigurator = fEmManager->EmConfigurator(); 144 emConfigurator->PrepareModels(&part, ptr, this); 145 } 146 147 if (fFirstParticle == &part) { 148 G4bool master = fEmManager->IsMaster(); 149 G4LossTableBuilder* bld = fEmManager->GetTableBuilder(); 150 G4bool baseMat = bld->GetBaseMaterialFlag(); 151 const auto* theParameters = G4EmParameters::Instance(); 152 153 if (master) { 154 SetVerboseLevel(theParameters->Verbose()); 155 } 156 else { 157 SetVerboseLevel(theParameters->WorkerVerbose()); 158 } 159 160 const G4int numberOfModels = fModelManager->NumberOfModels(); 161 if (fType == ScatteringType::MultipleScattering) { 162 for (G4int i = 0; i < numberOfModels; ++i) { 163 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i)); 164 msc->SetPolarAngleLimit(theParameters->MscThetaLimit()); 165 G4double emax = std::min(msc->HighEnergyLimit(), theParameters->MaxKinEnergy()); 166 msc->SetHighEnergyLimit(emax); 167 msc->SetUseBaseMaterials(baseMat); 168 } 169 } 170 else if (fType == ScatteringType::SingleScattering) { 171 if (master) { 172 if (fEmData == nullptr) { 173 fEmData = new G4EmDataHandler(2); 174 } 175 176 fLambdaTable = fEmData->MakeTable(0); 177 bld->InitialiseBaseMaterials(fLambdaTable); 178 } 179 } 180 181 fCuts = fModelManager->Initialise(fFirstParticle, G4Electron::Electron(), verboseLevel); 182 } 183 } 184 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 186 187 void G4TransportationWithMsc::BuildPhysicsTable(const G4ParticleDefinition& part) 188 { 189 if (fFirstParticle == &part) { 190 fEmManager->BuildPhysicsTable(fFirstParticle); 191 192 if (fEmManager->IsMaster()) { 193 if (fType == ScatteringType::SingleScattering) { 194 const auto* theParameters = G4EmParameters::Instance(); 195 G4LossTableBuilder* bld = fEmManager->GetTableBuilder(); 196 const G4ProductionCutsTable* theCoupleTable = 197 G4ProductionCutsTable::GetProductionCutsTable(); 198 std::size_t numOfCouples = theCoupleTable->GetTableSize(); 199 200 G4double emin = theParameters->MinKinEnergy(); 201 G4double emax = theParameters->MaxKinEnergy(); 202 203 G4double scale = emax / emin; 204 G4int nbin = theParameters->NumberOfBinsPerDecade() * G4lrint(std::log10(scale)); 205 scale = nbin / G4Log(scale); 206 207 G4int bin = G4lrint(scale * G4Log(emax / emin)); 208 bin = std::max(bin, 5); 209 210 for (std::size_t i = 0; i < numOfCouples; ++i) { 211 if (!bld->GetFlag(i)) continue; 212 213 // Create physics vector and fill it 214 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple((G4int)i); 215 216 auto* aVector = new G4PhysicsLogVector(emin, emax, bin, /*splineFlag*/ true); 217 fModelManager->FillLambdaVector(aVector, couple, /*startNull*/ false); 218 aVector->FillSecondDerivatives(); 219 G4PhysicsTableHelper::SetPhysicsVector(fLambdaTable, i, aVector); 220 } 221 } 222 } 223 else { 224 const auto masterProcess = static_cast<const G4TransportationWithMsc*>(GetMasterProcess()); 225 226 // Initialisation of models. 227 const G4int numberOfModels = fModelManager->NumberOfModels(); 228 if (fType == ScatteringType::MultipleScattering) { 229 for (G4int i = 0; i < numberOfModels; ++i) { 230 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i)); 231 auto msc0 = static_cast<G4VMscModel*>(masterProcess->fModelManager->GetModel(i)); 232 msc->SetCrossSectionTable(msc0->GetCrossSectionTable(), false); 233 msc->InitialiseLocal(fFirstParticle, msc0); 234 } 235 } 236 else if (fType == ScatteringType::SingleScattering) { 237 this->fLambdaTable = masterProcess->fLambdaTable; 238 } 239 } 240 } 241 242 if (!G4EmParameters::Instance()->IsPrintLocked() && verboseLevel > 0) { 243 G4cout << G4endl; 244 G4cout << GetProcessName() << ": for " << part.GetParticleName(); 245 if (fMultipleSteps) { 246 G4cout << " (multipleSteps: 1)"; 247 } 248 G4cout << G4endl; 249 fModelManager->DumpModelList(G4cout, verboseLevel); 250 } 251 } 252 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 254 255 void G4TransportationWithMsc::StartTracking(G4Track* track) 256 { 257 auto* currParticle = track->GetParticleDefinition(); 258 fIonisation = fEmManager->GetEnergyLossProcess(currParticle); 259 260 fSubStepDynamicParticle->SetDefinition(currParticle); 261 262 const G4int numberOfModels = fModelManager->NumberOfModels(); 263 if (fType == ScatteringType::MultipleScattering) { 264 for (G4int i = 0; i < numberOfModels; ++i) { 265 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i)); 266 msc->StartTracking(track); 267 msc->SetIonisation(fIonisation, currParticle); 268 } 269 } 270 271 // Ensure that field propagation state is also cleared / prepared 272 G4Transportation::StartTracking(track); 273 } 274 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 276 277 G4double G4TransportationWithMsc::AlongStepGetPhysicalInteractionLength(const G4Track& track, 278 G4double previousStepSize, 279 G4double currentMinimumStep, 280 G4double& proposedSafety, 281 G4GPILSelection* selection) 282 { 283 *selection = NotCandidateForSelection; 284 285 const G4double physStepLimit = currentMinimumStep; 286 287 switch (fType) { 288 case ScatteringType::MultipleScattering: { 289 // Select the MSC model for the current kinetic energy. 290 G4VMscModel* mscModel = nullptr; 291 const G4double ekin = track.GetKineticEnergy(); 292 const auto* couple = track.GetMaterialCutsCouple(); 293 const auto* particleDefinition = track.GetParticleDefinition(); 294 if (physStepLimit > kGeomMin) { 295 G4double ekinForSelection = ekin; 296 G4double pdgMass = particleDefinition->GetPDGMass(); 297 if (pdgMass > CLHEP::GeV) { 298 ekinForSelection *= proton_mass_c2 / pdgMass; 299 } 300 301 if (ekinForSelection >= kLowestKinEnergy) { 302 mscModel = static_cast<G4VMscModel*>( 303 fModelManager->SelectModel(ekinForSelection, couple->GetIndex())); 304 if (mscModel == nullptr) { 305 G4Exception("G4TransportationWithMsc::AlongStepGPIL", "em0052", FatalException, 306 "no MSC model found"); 307 } 308 if (!mscModel->IsActive(ekinForSelection)) { 309 mscModel = nullptr; 310 } 311 } 312 } 313 314 // Call the MSC model to potentially limit the step and convert to 315 // geometric path length. 316 if (mscModel != nullptr) { 317 mscModel->SetCurrentCouple(couple); 318 319 // Use the provided track for the first step. 320 const G4Track* currentTrackPtr = &track; 321 322 G4double currentSafety = proposedSafety; 323 G4double currentEnergy = ekin; 324 325 G4double stepLimitLeft = physStepLimit; 326 G4double totalGeometryStepLength = 0, totalTruePathLength = 0; 327 G4bool firstStep = true, continueStepping = fMultipleSteps; 328 329 do { 330 G4double gPathLength = stepLimitLeft; 331 G4double tPathLength = 332 mscModel->ComputeTruePathLengthLimit(*currentTrackPtr, gPathLength); 333 G4bool mscLimitsStep = (tPathLength < stepLimitLeft); 334 if (!fMultipleSteps && mscLimitsStep) { 335 // MSC limits the step. 336 *selection = CandidateForSelection; 337 } 338 339 if (!firstStep) { 340 // Move the navigator to where the previous step ended. 341 fLinearNavigator->LocateGlobalPointWithinVolume(fTransportEndPosition); 342 } 343 344 G4GPILSelection transportSelection; 345 G4double geometryStepLength = G4Transportation::AlongStepGetPhysicalInteractionLength( 346 *currentTrackPtr, previousStepSize, gPathLength, currentSafety, &transportSelection); 347 if (geometryStepLength < gPathLength) { 348 // Transportation limits the step, ie the track hit a boundary. 349 *selection = CandidateForSelection; 350 continueStepping = false; 351 } 352 if (fTransportEndKineticEnergy != currentEnergy) { 353 // Field propagation changed the energy, it's not possible to 354 // estimate the continuous energy loss and continue stepping. 355 continueStepping = false; 356 } 357 358 if (firstStep) { 359 proposedSafety = currentSafety; 360 } 361 totalGeometryStepLength += geometryStepLength; 362 363 // Sample MSC direction change and displacement. 364 const G4double range = mscModel->GetRange(particleDefinition, currentEnergy, couple); 365 366 tPathLength = mscModel->ComputeTrueStepLength(geometryStepLength); 367 368 // Protect against wrong t->g->t conversion. 369 tPathLength = std::min(tPathLength, stepLimitLeft); 370 371 totalTruePathLength += tPathLength; 372 if (*selection != CandidateForSelection && !mscLimitsStep) { 373 // If neither MSC nor transportation limits the step, we got the 374 // distance we want - make sure we exit the loop. 375 continueStepping = false; 376 } 377 else if (tPathLength >= range) { 378 // The particle will stop, exit the loop. 379 continueStepping = false; 380 } 381 else { 382 stepLimitLeft -= tPathLength; 383 } 384 385 // Do not sample scattering at the last or at a small step. 386 if (tPathLength < range && tPathLength > kGeomMin) { 387 static constexpr G4double minSafety = 1.20 * CLHEP::nm; 388 static constexpr G4double sFact = 0.99; 389 390 // The call to SampleScattering() *may* directly fill in the changed 391 // direction into fParticleChangeForMSC, so we have to: 392 // 1) Make sure the momentum direction is initialized. 393 fParticleChangeForMSC->ProposeMomentumDirection(fTransportEndMomentumDir); 394 // 2) Call SampleScattering(), which *may* change it. 395 const G4ThreeVector displacement = 396 mscModel->SampleScattering(fTransportEndMomentumDir, minSafety); 397 // 3) Get the changed direction. 398 fTransportEndMomentumDir = *fParticleChangeForMSC->GetProposedMomentumDirection(); 399 400 const G4double r2 = displacement.mag2(); 401 if (r2 > kMinDisplacement2) { 402 G4bool positionChanged = true; 403 G4double dispR = std::sqrt(r2); 404 G4double postSafety = 405 sFact * fpSafetyHelper->ComputeSafety(fTransportEndPosition, dispR); 406 407 // Far away from geometry boundary 408 if (postSafety > 0.0 && dispR <= postSafety) { 409 fTransportEndPosition += displacement; 410 411 // Near the boundary 412 } 413 else { 414 // displaced point is definitely within the volume 415 if (dispR < postSafety) { 416 fTransportEndPosition += displacement; 417 418 // reduced displacement 419 } 420 else if (postSafety > kGeomMin) { 421 fTransportEndPosition += displacement * (postSafety / dispR); 422 423 // very small postSafety 424 } 425 else { 426 positionChanged = false; 427 } 428 } 429 if (positionChanged) { 430 fpSafetyHelper->ReLocateWithinVolume(fTransportEndPosition); 431 } 432 } 433 } 434 435 if (continueStepping) { 436 // Update safety according to the geometry distance. 437 if (currentSafety < fEndPointDistance) { 438 currentSafety = 0; 439 } 440 else { 441 currentSafety -= fEndPointDistance; 442 } 443 444 // Update the kinetic energy according to the continuous loss. 445 currentEnergy = mscModel->GetEnergy(particleDefinition, range - tPathLength, couple); 446 447 // From now on, use the track that we can update below. 448 currentTrackPtr = fSubStepTrack; 449 450 fSubStepDynamicParticle->SetKineticEnergy(currentEnergy); 451 fSubStepDynamicParticle->SetMomentumDirection(fTransportEndMomentumDir); 452 fSubStepTrack->SetPosition(fTransportEndPosition); 453 454 G4StepPoint& subPreStepPoint = *fSubStep->GetPreStepPoint(); 455 subPreStepPoint.SetMaterialCutsCouple(couple); 456 subPreStepPoint.SetPosition(fTransportEndPosition); 457 subPreStepPoint.SetSafety(currentSafety); 458 subPreStepPoint.SetStepStatus(fAlongStepDoItProc); 459 } 460 firstStep = false; 461 } while (continueStepping); 462 463 // Note: currentEnergy is only updated if continueStepping is true. 464 // In case field propagation changed the energy, this flag is 465 // immediately set to false and currentEnergy is still equal to the 466 // initial kinetic energy stored in ekin. 467 if (currentEnergy != ekin) { 468 // If field propagation didn't change the energy and we potentially 469 // did multiple steps, reset the energy that G4Transportation will 470 // propose to not subtract the energy loss twice. 471 fTransportEndKineticEnergy = ekin; 472 // Also ask for the range again with the initial energy so it is 473 // correctly cached in the G4VEnergyLossProcess. 474 // FIXME: Asking for a range should never change the cached values! 475 (void)mscModel->GetRange(particleDefinition, ekin, couple); 476 } 477 478 fParticleChange.ProposeTrueStepLength(totalTruePathLength); 479 480 // Inform G4Transportation that the momentum might have changed due 481 // to scattering. We do this unconditionally to avoid the situation 482 // where the last step is done without MSC and G4Transportation reset 483 // the flag, for example when running without field. 484 fMomentumChanged = true; 485 486 return totalGeometryStepLength; 487 } 488 break; 489 } 490 491 case ScatteringType::SingleScattering: { 492 // Select the model for the current kinetic energy. 493 const G4double ekin = track.GetKineticEnergy(); 494 const auto* couple = track.GetMaterialCutsCouple(); 495 const auto* particleDefinition = track.GetParticleDefinition(); 496 497 G4double ekinForSelection = ekin; 498 G4double pdgMass = particleDefinition->GetPDGMass(); 499 if (pdgMass > CLHEP::GeV) { 500 ekinForSelection *= proton_mass_c2 / pdgMass; 501 } 502 503 G4VEmModel* currentModel = fModelManager->SelectModel(ekinForSelection, couple->GetIndex()); 504 if (currentModel == nullptr) { 505 G4Exception("G4TransportationWithMsc::AlongStepGPIL", "em0052", FatalException, 506 "no scattering model found"); 507 } 508 if (!currentModel->IsActive(ekinForSelection)) { 509 currentModel = nullptr; 510 } 511 512 if (currentModel != nullptr) { 513 currentModel->SetCurrentCouple(couple); 514 G4int coupleIndex = couple->GetIndex(); 515 516 // Compute mean free path. 517 G4double logEkin = track.GetDynamicParticle()->GetLogKineticEnergy(); 518 G4double lambda = ((*fLambdaTable)[coupleIndex])->LogVectorValue(ekin, logEkin); 519 if (lambda > 0.0) { 520 // Assume that the mean free path and dE/dx are constant along the 521 // step, which is a valid approximation for most cases. 522 G4double meanFreePath = 1.0 / lambda; 523 G4double dedx = fIonisation->GetDEDX(ekin, couple); 524 525 G4double currentSafety = proposedSafety; 526 G4double currentEnergy = ekin; 527 528 // Use the provided track for the first step. 529 const G4Track* currentTrackPtr = &track; 530 531 G4double stepLimitLeft = physStepLimit; 532 G4double totalStepLength = 0; 533 G4bool firstStep = true, continueStepping = fMultipleSteps; 534 535 do { 536 G4double interactionLength = meanFreePath * -G4Log(G4UniformRand()); 537 538 G4bool ssLimitsStep = (interactionLength < stepLimitLeft); 539 G4double gPathLength = stepLimitLeft; 540 if (ssLimitsStep) { 541 if (!fMultipleSteps) { 542 // Scattering limits the step. 543 *selection = CandidateForSelection; 544 } 545 gPathLength = interactionLength; 546 } 547 548 if (!firstStep) { 549 // Move the navigator to where the previous step ended. 550 fLinearNavigator->LocateGlobalPointWithinVolume(fTransportEndPosition); 551 } 552 553 G4GPILSelection transportSelection; 554 G4double geometryStepLength = G4Transportation::AlongStepGetPhysicalInteractionLength( 555 *currentTrackPtr, previousStepSize, gPathLength, currentSafety, &transportSelection); 556 if (geometryStepLength < gPathLength) { 557 // Transportation limits the step, ie the track hit a boundary. 558 *selection = CandidateForSelection; 559 ssLimitsStep = false; 560 continueStepping = false; 561 } 562 if (fTransportEndKineticEnergy != currentEnergy) { 563 // Field propagation changed the energy, it's not possible to 564 // estimate the continuous energy loss and continue stepping. 565 continueStepping = false; 566 } 567 568 if (firstStep) { 569 proposedSafety = currentSafety; 570 } 571 totalStepLength += geometryStepLength; 572 573 if (*selection != CandidateForSelection && !ssLimitsStep) { 574 // If neither scattering nor transportation limits the step, we 575 // got the distance we want - make sure we exit the loop. 576 continueStepping = false; 577 } 578 else { 579 stepLimitLeft -= geometryStepLength; 580 } 581 582 // Update the kinetic energy according to the continuous loss. 583 G4double energyAfterLinearLoss = 584 fTransportEndKineticEnergy - geometryStepLength * dedx; 585 586 if (ssLimitsStep) { 587 fSubStepDynamicParticle->SetKineticEnergy(energyAfterLinearLoss); 588 589 // The call to SampleSecondaries() directly fills in the changed 590 // direction into fParticleChangeForSS, so we have to: 591 // 1) Set the momentum direction in dynamic particle. 592 fSubStepDynamicParticle->SetMomentumDirection(fTransportEndMomentumDir); 593 // 2) Call SampleSecondaries(), which changes the direction. 594 currentModel->SampleSecondaries(fSecondariesSS, couple, fSubStepDynamicParticle, 595 (*fCuts)[coupleIndex]); 596 // 3) Get the changed direction. 597 fTransportEndMomentumDir = fParticleChangeForSS->GetProposedMomentumDirection(); 598 599 // Check that the model neither created secondaries nor proposed 600 // a local energy deposit because this process does not know how 601 // to handle these cases. 602 if (fSecondariesSS->size() > 0) { 603 G4Exception("G4TransportationWithMsc::AlongStepGPIL", "em0053", FatalException, 604 "scattering model created secondaries"); 605 } 606 if (fParticleChangeForSS->GetLocalEnergyDeposit() > 0) { 607 G4Exception("G4TransportationWithMsc::AlongStepGPIL", "em0053", FatalException, 608 "scattering model proposed energy deposit"); 609 } 610 } 611 612 if (continueStepping) { 613 // Update safety according to the geometry distance. 614 if (currentSafety < fEndPointDistance) { 615 currentSafety = 0; 616 } 617 else { 618 currentSafety -= fEndPointDistance; 619 } 620 621 // Update the energy taking continuous loss into account. 622 currentEnergy = energyAfterLinearLoss; 623 624 // From now on, use the track that we can update below. 625 currentTrackPtr = fSubStepTrack; 626 627 fSubStepDynamicParticle->SetKineticEnergy(currentEnergy); 628 fSubStepDynamicParticle->SetMomentumDirection(fTransportEndMomentumDir); 629 fSubStepTrack->SetPosition(fTransportEndPosition); 630 631 G4StepPoint& subPreStepPoint = *fSubStep->GetPreStepPoint(); 632 subPreStepPoint.SetMaterialCutsCouple(couple); 633 subPreStepPoint.SetPosition(fTransportEndPosition); 634 subPreStepPoint.SetSafety(currentSafety); 635 subPreStepPoint.SetStepStatus(fAlongStepDoItProc); 636 } 637 firstStep = false; 638 } while (continueStepping); 639 640 // Note: currentEnergy is only updated if continueStepping is true. 641 // In case field propagation changed the energy, this flag is 642 // immediately set to false and currentEnergy is still equal to the 643 // initial kinetic energy stored in ekin. 644 if (currentEnergy != ekin) { 645 // If field propagation didn't change the energy and we potentially 646 // did multiple steps, reset the energy that G4Transportation will 647 // propose to not subtract the energy loss twice. 648 fTransportEndKineticEnergy = ekin; 649 } 650 651 fParticleChange.ProposeTrueStepLength(totalStepLength); 652 653 // Inform G4Transportation that the momentum might have changed due 654 // to scattering, even if there is no field. 655 fMomentumChanged = true; 656 657 return totalStepLength; 658 } 659 } 660 break; 661 } 662 } 663 664 // If we get here, no scattering has happened. 665 return G4Transportation::AlongStepGetPhysicalInteractionLength( 666 track, previousStepSize, currentMinimumStep, proposedSafety, selection); 667 } 668 669 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 670