| 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"
>> 45 #include "G4VMscModel.hh"
>> 46
49 #include "G4ParticleChangeForMSC.hh" 47 #include "G4ParticleChangeForMSC.hh"
50 #include "G4PhysicalConstants.hh" << 48
51 #include "G4PhysicsTableHelper.hh" << 49 #include "G4DynamicParticle.hh"
52 #include "G4PhysicsVector.hh" <<
53 #include "G4ProductionCutsTable.hh" <<
54 #include "G4Step.hh" 50 #include "G4Step.hh"
55 #include "G4StepPoint.hh" 51 #include "G4StepPoint.hh"
56 #include "G4StepStatus.hh" 52 #include "G4StepStatus.hh"
57 #include "G4Track.hh" 53 #include "G4Track.hh"
58 #include "G4TransportationProcessType.hh" << 54
59 #include "G4VMscModel.hh" << 55 #include "G4Electron.hh"
>> 56
>> 57 #include "G4PhysicalConstants.hh"
60 58
61 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 60
63 static constexpr G4double kLowestKinEnergy = 1 << 61 static constexpr G4double kLowestKinEnergy = 10 * CLHEP::eV;
64 static constexpr G4double kGeomMin = 0.05 * CL << 62 static constexpr G4double kGeomMin = 0.05 * CLHEP::nm;
65 static constexpr G4double kMinDisplacement2 = 63 static constexpr G4double kMinDisplacement2 = kGeomMin * kGeomMin;
66 64
67 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 66
69 G4TransportationWithMsc::G4TransportationWithM << 67 G4TransportationWithMsc::G4TransportationWithMsc(ScatteringType type,
70 : G4Transportation(verbosity, "Transportatio << 68 G4int verbosity)
>> 69 : G4Transportation(verbosity, "TransportationWithMsc")
>> 70 , fType(type)
71 { 71 {
72 SetProcessSubType(static_cast<G4int>(TRANSPO <<
73 SetVerboseLevel(1); 72 SetVerboseLevel(1);
74 73
75 fEmManager = G4LossTableManager::Instance(); << 74 fEmManager = G4LossTableManager::Instance();
76 fModelManager = new G4EmModelManager; 75 fModelManager = new G4EmModelManager;
77 76
78 if (type == ScatteringType::MultipleScatteri << 77 if(type == ScatteringType::MultipleScattering)
>> 78 {
79 fParticleChangeForMSC = new G4ParticleChan 79 fParticleChangeForMSC = new G4ParticleChangeForMSC;
80 } 80 }
81 else if (type == ScatteringType::SingleScatt <<
82 fParticleChangeForSS = new G4ParticleChang <<
83 fSecondariesSS = new std::vector<G4Dynamic <<
84 } <<
85 81
86 G4ThreeVector zero; 82 G4ThreeVector zero;
87 fSubStepDynamicParticle = new G4DynamicParti << 83 fSubStepDynamicParticle =
>> 84 new G4DynamicParticle(G4Electron::Definition(), zero);
88 fSubStepTrack = new G4Track(fSubStepDynamicP 85 fSubStepTrack = new G4Track(fSubStepDynamicParticle, 0, zero);
89 fSubStep = new G4Step; << 86 fSubStep = new G4Step;
90 fSubStepTrack->SetStep(fSubStep); 87 fSubStepTrack->SetStep(fSubStep);
91 } 88 }
92 89
93 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 91
95 G4TransportationWithMsc::~G4TransportationWith 92 G4TransportationWithMsc::~G4TransportationWithMsc()
96 { 93 {
97 delete fModelManager; 94 delete fModelManager;
98 delete fParticleChangeForMSC; 95 delete fParticleChangeForMSC;
99 delete fEmData; <<
100 delete fParticleChangeForSS; <<
101 delete fSecondariesSS; <<
102 96
103 // fSubStepDynamicParticle is owned and also 97 // fSubStepDynamicParticle is owned and also deleted by fSubStepTrack!
104 delete fSubStepTrack; 98 delete fSubStepTrack;
105 delete fSubStep; 99 delete fSubStep;
106 } 100 }
107 101
108 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 103
110 void G4TransportationWithMsc::AddMscModel(G4VM 104 void G4TransportationWithMsc::AddMscModel(G4VMscModel* mscModel, G4int order,
111 cons 105 const G4Region* region)
112 { 106 {
113 if (fType != ScatteringType::MultipleScatter << 107 if(fType != ScatteringType::MultipleScattering)
114 G4Exception("G4TransportationWithMsc::AddM << 108 {
>> 109 G4Exception("G4TransportationWithMsc::AddMscModel", "em0051",
>> 110 FatalException,
115 "not allowed unless type == Mu 111 "not allowed unless type == MultipleScattering");
116 } 112 }
117 113
118 fModelManager->AddEmModel(order, mscModel, n 114 fModelManager->AddEmModel(order, mscModel, nullptr, region);
119 mscModel->SetParticleChange(fParticleChangeF 115 mscModel->SetParticleChange(fParticleChangeForMSC);
120 } 116 }
121 117
122 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 119
124 void G4TransportationWithMsc::AddSSModel(G4VEm << 120 void G4TransportationWithMsc::PreparePhysicsTable(
>> 121 const G4ParticleDefinition& part)
125 { 122 {
126 if (fType != ScatteringType::SingleScatterin << 123 if(nullptr == fFirstParticle)
127 G4Exception("G4TransportationWithMsc::AddS << 124 {
128 "not allowed unless type == Si <<
129 } <<
130 <<
131 fModelManager->AddEmModel(order, model, null <<
132 model->SetPolarAngleLimit(0.0); <<
133 model->SetParticleChange(fParticleChangeForS <<
134 } <<
135 <<
136 //....oooOO0OOooo........oooOO0OOooo........oo <<
137 <<
138 void G4TransportationWithMsc::PreparePhysicsTa <<
139 { <<
140 if (nullptr == fFirstParticle) { <<
141 fFirstParticle = ∂ 125 fFirstParticle = ∂
142 G4VMultipleScattering* ptr = nullptr; 126 G4VMultipleScattering* ptr = nullptr;
143 auto emConfigurator = fEmManager->EmConfig 127 auto emConfigurator = fEmManager->EmConfigurator();
144 emConfigurator->PrepareModels(&part, ptr, 128 emConfigurator->PrepareModels(&part, ptr, this);
145 } 129 }
146 130
147 if (fFirstParticle == &part) { << 131 if(fFirstParticle == &part)
148 G4bool master = fEmManager->IsMaster(); << 132 {
149 G4LossTableBuilder* bld = fEmManager->GetT << 133 G4bool master = fEmManager->IsMaster();
150 G4bool baseMat = bld->GetBaseMaterialFlag( << 134 G4LossTableBuilder* bld = fEmManager->GetTableBuilder();
>> 135 G4bool baseMat = bld->GetBaseMaterialFlag();
151 const auto* theParameters = G4EmParameters 136 const auto* theParameters = G4EmParameters::Instance();
152 137
153 if (master) { << 138 if(master)
>> 139 {
154 SetVerboseLevel(theParameters->Verbose() 140 SetVerboseLevel(theParameters->Verbose());
155 } 141 }
156 else { << 142 else
>> 143 {
157 SetVerboseLevel(theParameters->WorkerVer 144 SetVerboseLevel(theParameters->WorkerVerbose());
158 } 145 }
159 146
160 const G4int numberOfModels = fModelManager 147 const G4int numberOfModels = fModelManager->NumberOfModels();
161 if (fType == ScatteringType::MultipleScatt << 148 if(fType == ScatteringType::MultipleScattering)
162 for (G4int i = 0; i < numberOfModels; ++ << 149 {
>> 150 for(G4int i = 0; i < numberOfModels; ++i)
>> 151 {
163 auto msc = static_cast<G4VMscModel*>(f 152 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i));
>> 153 msc->SetMasterThread(master);
164 msc->SetPolarAngleLimit(theParameters- 154 msc->SetPolarAngleLimit(theParameters->MscThetaLimit());
165 G4double emax = std::min(msc->HighEner << 155 G4double emax =
>> 156 std::min(msc->HighEnergyLimit(), theParameters->MaxKinEnergy());
166 msc->SetHighEnergyLimit(emax); 157 msc->SetHighEnergyLimit(emax);
167 msc->SetUseBaseMaterials(baseMat); 158 msc->SetUseBaseMaterials(baseMat);
168 } 159 }
169 } 160 }
170 else if (fType == ScatteringType::SingleSc <<
171 if (master) { <<
172 if (fEmData == nullptr) { <<
173 fEmData = new G4EmDataHandler(2); <<
174 } <<
175 161
176 fLambdaTable = fEmData->MakeTable(0); << 162 fModelManager->Initialise(fFirstParticle, G4Electron::Electron(),
177 bld->InitialiseBaseMaterials(fLambdaTa << 163 verboseLevel);
178 } <<
179 } <<
180 <<
181 fCuts = fModelManager->Initialise(fFirstPa <<
182 } 164 }
183 } 165 }
184 166
185 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186 168
187 void G4TransportationWithMsc::BuildPhysicsTabl << 169 void G4TransportationWithMsc::BuildPhysicsTable(
>> 170 const G4ParticleDefinition& part)
188 { 171 {
189 if (fFirstParticle == &part) { << 172 if(fFirstParticle == &part)
>> 173 {
190 fEmManager->BuildPhysicsTable(fFirstPartic 174 fEmManager->BuildPhysicsTable(fFirstParticle);
191 175
192 if (fEmManager->IsMaster()) { << 176 if(!fEmManager->IsMaster())
193 if (fType == ScatteringType::SingleScatt << 177 {
194 const auto* theParameters = G4EmParame << 178 const auto masterProcess =
195 G4LossTableBuilder* bld = fEmManager-> << 179 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 180
226 // Initialisation of models. 181 // Initialisation of models.
227 const G4int numberOfModels = fModelManag 182 const G4int numberOfModels = fModelManager->NumberOfModels();
228 if (fType == ScatteringType::MultipleSca << 183 if(fType == ScatteringType::MultipleScattering)
229 for (G4int i = 0; i < numberOfModels; << 184 {
>> 185 for(G4int i = 0; i < numberOfModels; ++i)
>> 186 {
230 auto msc = static_cast<G4VMscModel*> 187 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i));
231 auto msc0 = static_cast<G4VMscModel* << 188 auto msc0 =
>> 189 static_cast<G4VMscModel*>(masterProcess->fModelManager->GetModel(i));
232 msc->SetCrossSectionTable(msc0->GetC 190 msc->SetCrossSectionTable(msc0->GetCrossSectionTable(), false);
233 msc->InitialiseLocal(fFirstParticle, 191 msc->InitialiseLocal(fFirstParticle, msc0);
234 } 192 }
235 } 193 }
236 else if (fType == ScatteringType::Single <<
237 this->fLambdaTable = masterProcess->fL <<
238 } <<
239 } 194 }
240 } 195 }
241 196
242 if (!G4EmParameters::Instance()->IsPrintLock << 197 if(!G4EmParameters::Instance()->IsPrintLocked() && verboseLevel > 0)
>> 198 {
243 G4cout << G4endl; 199 G4cout << G4endl;
244 G4cout << GetProcessName() << ": for " << 200 G4cout << GetProcessName() << ": for " << part.GetParticleName();
245 if (fMultipleSteps) { << 201 if(fMultipleSteps)
>> 202 {
246 G4cout << " (multipleSteps: 1)"; 203 G4cout << " (multipleSteps: 1)";
247 } 204 }
248 G4cout << G4endl; 205 G4cout << G4endl;
249 fModelManager->DumpModelList(G4cout, verbo 206 fModelManager->DumpModelList(G4cout, verboseLevel);
250 } 207 }
251 } 208 }
252 209
253 //....oooOO0OOooo........oooOO0OOooo........oo 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 211
255 void G4TransportationWithMsc::StartTracking(G4 212 void G4TransportationWithMsc::StartTracking(G4Track* track)
256 { 213 {
257 auto* currParticle = track->GetParticleDefin 214 auto* currParticle = track->GetParticleDefinition();
258 fIonisation = fEmManager->GetEnergyLossProce << 215 auto* ionisation = fEmManager->GetEnergyLossProcess(currParticle);
259 216
260 fSubStepDynamicParticle->SetDefinition(currP 217 fSubStepDynamicParticle->SetDefinition(currParticle);
261 218
262 const G4int numberOfModels = fModelManager-> 219 const G4int numberOfModels = fModelManager->NumberOfModels();
263 if (fType == ScatteringType::MultipleScatter << 220 if(fType == ScatteringType::MultipleScattering)
264 for (G4int i = 0; i < numberOfModels; ++i) << 221 {
>> 222 for(G4int i = 0; i < numberOfModels; ++i)
>> 223 {
265 auto msc = static_cast<G4VMscModel*>(fMo 224 auto msc = static_cast<G4VMscModel*>(fModelManager->GetModel(i));
266 msc->StartTracking(track); 225 msc->StartTracking(track);
267 msc->SetIonisation(fIonisation, currPart << 226 msc->SetIonisation(ionisation, currParticle);
268 } 227 }
269 } 228 }
270 229
271 // Ensure that field propagation state is al 230 // Ensure that field propagation state is also cleared / prepared
272 G4Transportation::StartTracking(track); 231 G4Transportation::StartTracking(track);
273 } 232 }
274 233
275 //....oooOO0OOooo........oooOO0OOooo........oo 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 235
277 G4double G4TransportationWithMsc::AlongStepGet << 236 G4double G4TransportationWithMsc::AlongStepGetPhysicalInteractionLength(
278 << 237 const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
279 << 238 G4double& proposedSafety, G4GPILSelection* selection)
280 <<
281 <<
282 { 239 {
283 *selection = NotCandidateForSelection; 240 *selection = NotCandidateForSelection;
284 241
285 const G4double physStepLimit = currentMinimu 242 const G4double physStepLimit = currentMinimumStep;
286 243
287 switch (fType) { << 244 switch(fType)
>> 245 {
288 case ScatteringType::MultipleScattering: { 246 case ScatteringType::MultipleScattering: {
289 // Select the MSC model for the current 247 // Select the MSC model for the current kinetic energy.
290 G4VMscModel* mscModel = nullptr; << 248 G4VMscModel* mscModel = nullptr;
291 const G4double ekin = track.GetKineticEn << 249 const G4double ekin = track.GetKineticEnergy();
292 const auto* couple = track.GetMaterialCu << 250 const auto* couple = track.GetMaterialCutsCouple();
293 const auto* particleDefinition = track.G 251 const auto* particleDefinition = track.GetParticleDefinition();
294 if (physStepLimit > kGeomMin) { << 252 if(physStepLimit > kGeomMin)
>> 253 {
295 G4double ekinForSelection = ekin; 254 G4double ekinForSelection = ekin;
296 G4double pdgMass = particleDefinition- << 255 G4double pdgMass = particleDefinition->GetPDGMass();
297 if (pdgMass > CLHEP::GeV) { << 256 if(pdgMass > CLHEP::GeV)
>> 257 {
298 ekinForSelection *= proton_mass_c2 / 258 ekinForSelection *= proton_mass_c2 / pdgMass;
299 } 259 }
300 260
301 if (ekinForSelection >= kLowestKinEner << 261 if(ekinForSelection >= kLowestKinEnergy)
>> 262 {
302 mscModel = static_cast<G4VMscModel*> 263 mscModel = static_cast<G4VMscModel*>(
303 fModelManager->SelectModel(ekinFor 264 fModelManager->SelectModel(ekinForSelection, couple->GetIndex()));
304 if (mscModel == nullptr) { << 265 if(mscModel == nullptr)
305 G4Exception("G4TransportationWithM << 266 {
306 "no MSC model found"); << 267 G4Exception("G4TransportationWithMsc::AlongStepGPIL", "em0052",
>> 268 FatalException, "no MSC model found");
307 } 269 }
308 if (!mscModel->IsActive(ekinForSelec << 270 if(!mscModel->IsActive(ekinForSelection))
>> 271 {
309 mscModel = nullptr; 272 mscModel = nullptr;
310 } 273 }
311 } 274 }
312 } 275 }
313 276
314 // Call the MSC model to potentially lim 277 // Call the MSC model to potentially limit the step and convert to
315 // geometric path length. 278 // geometric path length.
316 if (mscModel != nullptr) { << 279 if(mscModel != nullptr)
>> 280 {
317 mscModel->SetCurrentCouple(couple); 281 mscModel->SetCurrentCouple(couple);
318 282
319 // Use the provided track for the firs 283 // Use the provided track for the first step.
320 const G4Track* currentTrackPtr = &trac 284 const G4Track* currentTrackPtr = &track;
321 285
322 G4double currentSafety = proposedSafet 286 G4double currentSafety = proposedSafety;
323 G4double currentEnergy = ekin; 287 G4double currentEnergy = ekin;
324 288
325 G4double stepLimitLeft = physStepLimit << 289 G4double stepLimitLeft = physStepLimit;
326 G4double totalGeometryStepLength = 0, 290 G4double totalGeometryStepLength = 0, totalTruePathLength = 0;
327 G4bool firstStep = true, continueStepp 291 G4bool firstStep = true, continueStepping = fMultipleSteps;
328 292
329 do { << 293 do
>> 294 {
330 G4double gPathLength = stepLimitLeft 295 G4double gPathLength = stepLimitLeft;
331 G4double tPathLength = 296 G4double tPathLength =
332 mscModel->ComputeTruePathLengthLim 297 mscModel->ComputeTruePathLengthLimit(*currentTrackPtr, gPathLength);
333 G4bool mscLimitsStep = (tPathLength 298 G4bool mscLimitsStep = (tPathLength < stepLimitLeft);
334 if (!fMultipleSteps && mscLimitsStep << 299 if(!fMultipleSteps && mscLimitsStep)
>> 300 {
335 // MSC limits the step. 301 // MSC limits the step.
336 *selection = CandidateForSelection 302 *selection = CandidateForSelection;
337 } 303 }
338 304
339 if (!firstStep) { << 305 if(!firstStep)
>> 306 {
340 // Move the navigator to where the 307 // Move the navigator to where the previous step ended.
341 fLinearNavigator->LocateGlobalPoin << 308 fLinearNavigator->LocateGlobalPointWithinVolume(
>> 309 fTransportEndPosition);
342 } 310 }
343 311
344 G4GPILSelection transportSelection; 312 G4GPILSelection transportSelection;
345 G4double geometryStepLength = G4Tran << 313 G4double geometryStepLength =
346 *currentTrackPtr, previousStepSize << 314 G4Transportation::AlongStepGetPhysicalInteractionLength(
347 if (geometryStepLength < gPathLength << 315 *currentTrackPtr, previousStepSize, gPathLength, currentSafety,
>> 316 &transportSelection);
>> 317 if(geometryStepLength < gPathLength)
>> 318 {
348 // Transportation limits the step, 319 // Transportation limits the step, ie the track hit a boundary.
349 *selection = CandidateForSelection << 320 *selection = CandidateForSelection;
350 continueStepping = false; 321 continueStepping = false;
351 } 322 }
352 if (fTransportEndKineticEnergy != cu << 323 if(fTransportEndKineticEnergy != currentEnergy)
>> 324 {
353 // Field propagation changed the e 325 // Field propagation changed the energy, it's not possible to
354 // estimate the continuous energy 326 // estimate the continuous energy loss and continue stepping.
355 continueStepping = false; 327 continueStepping = false;
356 } 328 }
357 329
358 if (firstStep) { << 330 if(firstStep)
>> 331 {
359 proposedSafety = currentSafety; 332 proposedSafety = currentSafety;
360 } 333 }
361 totalGeometryStepLength += geometryS 334 totalGeometryStepLength += geometryStepLength;
362 335
363 // Sample MSC direction change and d 336 // Sample MSC direction change and displacement.
364 const G4double range = mscModel->Get << 337 const G4double range =
>> 338 mscModel->GetRange(particleDefinition, currentEnergy, couple);
365 339
366 tPathLength = mscModel->ComputeTrueS 340 tPathLength = mscModel->ComputeTrueStepLength(geometryStepLength);
367 341
368 // Protect against wrong t->g->t con 342 // Protect against wrong t->g->t conversion.
369 tPathLength = std::min(tPathLength, 343 tPathLength = std::min(tPathLength, stepLimitLeft);
370 344
371 totalTruePathLength += tPathLength; 345 totalTruePathLength += tPathLength;
372 if (*selection != CandidateForSelect << 346 if(*selection != CandidateForSelection && !mscLimitsStep)
>> 347 {
373 // If neither MSC nor transportati 348 // If neither MSC nor transportation limits the step, we got the
374 // distance we want - make sure we 349 // distance we want - make sure we exit the loop.
375 continueStepping = false; 350 continueStepping = false;
376 } 351 }
377 else if (tPathLength >= range) { << 352 else if(tPathLength >= range)
>> 353 {
378 // The particle will stop, exit th 354 // The particle will stop, exit the loop.
379 continueStepping = false; 355 continueStepping = false;
380 } 356 }
381 else { << 357 else
>> 358 {
382 stepLimitLeft -= tPathLength; 359 stepLimitLeft -= tPathLength;
383 } 360 }
384 361
385 // Do not sample scattering at the l 362 // Do not sample scattering at the last or at a small step.
386 if (tPathLength < range && tPathLeng << 363 if(tPathLength < range && tPathLength > kGeomMin)
>> 364 {
387 static constexpr G4double minSafet 365 static constexpr G4double minSafety = 1.20 * CLHEP::nm;
388 static constexpr G4double sFact = << 366 static constexpr G4double sFact = 0.99;
389 367
390 // The call to SampleScattering() 368 // The call to SampleScattering() *may* directly fill in the changed
391 // direction into fParticleChangeF 369 // direction into fParticleChangeForMSC, so we have to:
392 // 1) Make sure the momentum direc 370 // 1) Make sure the momentum direction is initialized.
393 fParticleChangeForMSC->ProposeMome 371 fParticleChangeForMSC->ProposeMomentumDirection(fTransportEndMomentumDir);
394 // 2) Call SampleScattering(), whi 372 // 2) Call SampleScattering(), which *may* change it.
395 const G4ThreeVector displacement = 373 const G4ThreeVector displacement =
396 mscModel->SampleScattering(fTran 374 mscModel->SampleScattering(fTransportEndMomentumDir, minSafety);
397 // 3) Get the changed direction. 375 // 3) Get the changed direction.
398 fTransportEndMomentumDir = *fParti 376 fTransportEndMomentumDir = *fParticleChangeForMSC->GetProposedMomentumDirection();
399 377
400 const G4double r2 = displacement.m 378 const G4double r2 = displacement.mag2();
401 if (r2 > kMinDisplacement2) { << 379 if(r2 > kMinDisplacement2)
>> 380 {
402 G4bool positionChanged = true; 381 G4bool positionChanged = true;
403 G4double dispR = std::sqrt(r2); << 382 G4double dispR = std::sqrt(r2);
404 G4double postSafety = << 383 G4double postSafety = sFact * fpSafetyHelper->ComputeSafety(
405 sFact * fpSafetyHelper->Comput << 384 fTransportEndPosition, dispR);
406 385
407 // Far away from geometry bounda 386 // Far away from geometry boundary
408 if (postSafety > 0.0 && dispR <= << 387 if(postSafety > 0.0 && dispR <= postSafety)
>> 388 {
409 fTransportEndPosition += displ 389 fTransportEndPosition += displacement;
410 390
411 // Near the boundary 391 // Near the boundary
412 } 392 }
413 else { << 393 else
>> 394 {
414 // displaced point is definite 395 // displaced point is definitely within the volume
415 if (dispR < postSafety) { << 396 if(dispR < postSafety)
>> 397 {
416 fTransportEndPosition += dis 398 fTransportEndPosition += displacement;
417 399
418 // reduced displacement 400 // reduced displacement
419 } 401 }
420 else if (postSafety > kGeomMin << 402 else if(postSafety > kGeomMin)
>> 403 {
421 fTransportEndPosition += dis 404 fTransportEndPosition += displacement * (postSafety / dispR);
422 405
423 // very small postSafety 406 // very small postSafety
424 } 407 }
425 else { << 408 else
>> 409 {
426 positionChanged = false; 410 positionChanged = false;
427 } 411 }
428 } 412 }
429 if (positionChanged) { << 413 if(positionChanged)
>> 414 {
430 fpSafetyHelper->ReLocateWithin 415 fpSafetyHelper->ReLocateWithinVolume(fTransportEndPosition);
431 } 416 }
432 } 417 }
433 } 418 }
434 419
435 if (continueStepping) { << 420 if(continueStepping)
>> 421 {
436 // Update safety according to the 422 // Update safety according to the geometry distance.
437 if (currentSafety < fEndPointDista << 423 if(currentSafety < fEndPointDistance)
>> 424 {
438 currentSafety = 0; 425 currentSafety = 0;
439 } 426 }
440 else { << 427 else
>> 428 {
441 currentSafety -= fEndPointDistan 429 currentSafety -= fEndPointDistance;
442 } 430 }
443 431
444 // Update the kinetic energy accor 432 // Update the kinetic energy according to the continuous loss.
445 currentEnergy = mscModel->GetEnerg << 433 currentEnergy = mscModel->GetEnergy(particleDefinition,
>> 434 range - tPathLength, couple);
446 435
447 // From now on, use the track that 436 // From now on, use the track that we can update below.
448 currentTrackPtr = fSubStepTrack; 437 currentTrackPtr = fSubStepTrack;
449 438
450 fSubStepDynamicParticle->SetKineti 439 fSubStepDynamicParticle->SetKineticEnergy(currentEnergy);
451 fSubStepDynamicParticle->SetMoment << 440 fSubStepDynamicParticle->SetMomentumDirection(
>> 441 fTransportEndMomentumDir);
452 fSubStepTrack->SetPosition(fTransp 442 fSubStepTrack->SetPosition(fTransportEndPosition);
453 443
454 G4StepPoint& subPreStepPoint = *fS 444 G4StepPoint& subPreStepPoint = *fSubStep->GetPreStepPoint();
455 subPreStepPoint.SetMaterialCutsCou 445 subPreStepPoint.SetMaterialCutsCouple(couple);
456 subPreStepPoint.SetPosition(fTrans 446 subPreStepPoint.SetPosition(fTransportEndPosition);
457 subPreStepPoint.SetSafety(currentS 447 subPreStepPoint.SetSafety(currentSafety);
458 subPreStepPoint.SetStepStatus(fAlo 448 subPreStepPoint.SetStepStatus(fAlongStepDoItProc);
459 } 449 }
460 firstStep = false; 450 firstStep = false;
461 } while (continueStepping); << 451 } while(continueStepping);
462 452
463 // Note: currentEnergy is only updated 453 // Note: currentEnergy is only updated if continueStepping is true.
464 // In case field propagation changed t 454 // In case field propagation changed the energy, this flag is
465 // immediately set to false and curren 455 // immediately set to false and currentEnergy is still equal to the
466 // initial kinetic energy stored in ek 456 // initial kinetic energy stored in ekin.
467 if (currentEnergy != ekin) { << 457 if(currentEnergy != ekin)
>> 458 {
468 // If field propagation didn't chang 459 // If field propagation didn't change the energy and we potentially
469 // did multiple steps, reset the ene 460 // did multiple steps, reset the energy that G4Transportation will
470 // propose to not subtract the energ 461 // propose to not subtract the energy loss twice.
471 fTransportEndKineticEnergy = ekin; 462 fTransportEndKineticEnergy = ekin;
472 // Also ask for the range again with 463 // Also ask for the range again with the initial energy so it is
473 // correctly cached in the G4VEnergy 464 // correctly cached in the G4VEnergyLossProcess.
474 // FIXME: Asking for a range should 465 // FIXME: Asking for a range should never change the cached values!
475 (void)mscModel->GetRange(particleDef << 466 (void) mscModel->GetRange(particleDefinition, ekin, couple);
476 } 467 }
477 468
478 fParticleChange.ProposeTrueStepLength( 469 fParticleChange.ProposeTrueStepLength(totalTruePathLength);
479 470
480 // Inform G4Transportation that the mo 471 // Inform G4Transportation that the momentum might have changed due
481 // to scattering. We do this unconditi 472 // to scattering. We do this unconditionally to avoid the situation
482 // where the last step is done without 473 // where the last step is done without MSC and G4Transportation reset
483 // the flag, for example when running 474 // the flag, for example when running without field.
484 fMomentumChanged = true; 475 fMomentumChanged = true;
485 476
486 return totalGeometryStepLength; 477 return totalGeometryStepLength;
487 } 478 }
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 } 479 }
662 } 480 }
663 481
664 // If we get here, no scattering has happene 482 // If we get here, no scattering has happened.
665 return G4Transportation::AlongStepGetPhysica 483 return G4Transportation::AlongStepGetPhysicalInteractionLength(
666 track, previousStepSize, currentMinimumSte 484 track, previousStepSize, currentMinimumStep, proposedSafety, selection);
667 } 485 }
668 486
669 //....oooOO0OOooo........oooOO0OOooo........oo 487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
670 488