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