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