| 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 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // GEANT4 Class file 28 // GEANT4 Class file
29 // 29 //
30 // 30 //
31 // File name: G4VEnergyLossProcess 31 // File name: G4VEnergyLossProcess
32 // 32 //
33 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko
34 // 34 //
35 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002
36 // 36 //
37 // Modifications: Vladimir Ivanchenko 37 // Modifications: Vladimir Ivanchenko
38 // 38 //
39 // 39 //
40 // Class Description: 40 // Class Description:
41 // 41 //
42 // It is the unified energy loss process it ca 42 // It is the unified energy loss process it calculates the continuous
43 // energy loss for charged particles using a s 43 // energy loss for charged particles using a set of Energy Loss
44 // models valid for different energy regions. 44 // models valid for different energy regions. There are a possibility
45 // to create and access to dE/dx and range tab 45 // to create and access to dE/dx and range tables, or to calculate
46 // that information on fly. 46 // that information on fly.
47 // ------------------------------------------- 47 // -------------------------------------------------------------------
48 // 48 //
49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 51
52 #include "G4VEnergyLossProcess.hh" 52 #include "G4VEnergyLossProcess.hh"
53 #include "G4PhysicalConstants.hh" 53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh" 54 #include "G4SystemOfUnits.hh"
55 #include "G4ProcessManager.hh" 55 #include "G4ProcessManager.hh"
56 #include "G4LossTableManager.hh" 56 #include "G4LossTableManager.hh"
57 #include "G4LossTableBuilder.hh" 57 #include "G4LossTableBuilder.hh"
58 #include "G4Step.hh" 58 #include "G4Step.hh"
59 #include "G4ParticleDefinition.hh" 59 #include "G4ParticleDefinition.hh"
60 #include "G4ParticleTable.hh" 60 #include "G4ParticleTable.hh"
61 #include "G4EmParameters.hh" 61 #include "G4EmParameters.hh"
62 #include "G4EmUtility.hh" 62 #include "G4EmUtility.hh"
63 #include "G4EmTableUtil.hh" 63 #include "G4EmTableUtil.hh"
64 #include "G4VEmModel.hh" 64 #include "G4VEmModel.hh"
65 #include "G4VEmFluctuationModel.hh" 65 #include "G4VEmFluctuationModel.hh"
66 #include "G4DataVector.hh" 66 #include "G4DataVector.hh"
67 #include "G4PhysicsLogVector.hh" 67 #include "G4PhysicsLogVector.hh"
68 #include "G4VParticleChange.hh" 68 #include "G4VParticleChange.hh"
69 #include "G4Electron.hh" 69 #include "G4Electron.hh"
70 #include "G4ProcessManager.hh" 70 #include "G4ProcessManager.hh"
71 #include "G4UnitsTable.hh" 71 #include "G4UnitsTable.hh"
72 #include "G4Region.hh" 72 #include "G4Region.hh"
73 #include "G4RegionStore.hh" 73 #include "G4RegionStore.hh"
74 #include "G4PhysicsTableHelper.hh" 74 #include "G4PhysicsTableHelper.hh"
75 #include "G4SafetyHelper.hh" 75 #include "G4SafetyHelper.hh"
76 #include "G4EmDataHandler.hh" 76 #include "G4EmDataHandler.hh"
77 #include "G4TransportationManager.hh" 77 #include "G4TransportationManager.hh"
78 #include "G4VAtomDeexcitation.hh" 78 #include "G4VAtomDeexcitation.hh"
79 #include "G4VSubCutProducer.hh" 79 #include "G4VSubCutProducer.hh"
80 #include "G4EmBiasingManager.hh" 80 #include "G4EmBiasingManager.hh"
81 #include "G4Log.hh" 81 #include "G4Log.hh"
82 #include <iostream> 82 #include <iostream>
83 83
84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 85
86 namespace 86 namespace
87 { 87 {
88 G4String tnames[7] = 88 G4String tnames[7] =
89 {"DEDX","Ionisation","DEDXnr","CSDARange", 89 {"DEDX","Ionisation","DEDXnr","CSDARange","Lambda","Range","InverseRange"};
90 } 90 }
91 91
92 92
93 G4VEnergyLossProcess::G4VEnergyLossProcess(con 93 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name,
94 G4P 94 G4ProcessType type):
95 G4VContinuousDiscreteProcess(name, type) 95 G4VContinuousDiscreteProcess(name, type)
96 { 96 {
97 theParameters = G4EmParameters::Instance(); 97 theParameters = G4EmParameters::Instance();
98 SetVerboseLevel(1); 98 SetVerboseLevel(1);
99 99
100 // low energy limit 100 // low energy limit
101 lowestKinEnergy = theParameters->LowestElect 101 lowestKinEnergy = theParameters->LowestElectronEnergy();
102 102
103 // Size of tables 103 // Size of tables
104 minKinEnergy = 0.1*CLHEP::keV; 104 minKinEnergy = 0.1*CLHEP::keV;
105 maxKinEnergy = 100.0*CLHEP::TeV; 105 maxKinEnergy = 100.0*CLHEP::TeV;
106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; 106 maxKinEnergyCSDA = 1.0*CLHEP::GeV;
107 nBins = 84; 107 nBins = 84;
108 nBinsCSDA = 35; 108 nBinsCSDA = 35;
109 109
110 invLambdaFactor = 1.0/lambdaFactor; 110 invLambdaFactor = 1.0/lambdaFactor;
111 111
112 // default linear loss limit 112 // default linear loss limit
113 finalRange = 1.*CLHEP::mm; 113 finalRange = 1.*CLHEP::mm;
114 114
115 // run time objects 115 // run time objects
116 pParticleChange = &fParticleChange; 116 pParticleChange = &fParticleChange;
117 fParticleChange.SetSecondaryWeightByProcess( 117 fParticleChange.SetSecondaryWeightByProcess(true);
118 modelManager = new G4EmModelManager(); 118 modelManager = new G4EmModelManager();
119 safetyHelper = G4TransportationManager::GetT 119 safetyHelper = G4TransportationManager::GetTransportationManager()
120 ->GetSafetyHelper(); 120 ->GetSafetyHelper();
121 aGPILSelection = CandidateForSelection; 121 aGPILSelection = CandidateForSelection;
122 122
123 // initialise model 123 // initialise model
124 lManager = G4LossTableManager::Instance(); 124 lManager = G4LossTableManager::Instance();
125 lManager->Register(this); 125 lManager->Register(this);
126 isMaster = lManager->IsMaster(); 126 isMaster = lManager->IsMaster();
127 127
128 G4LossTableBuilder* bld = lManager->GetTable 128 G4LossTableBuilder* bld = lManager->GetTableBuilder();
129 theDensityFactor = bld->GetDensityFactors(); 129 theDensityFactor = bld->GetDensityFactors();
130 theDensityIdx = bld->GetCoupleIndexes(); 130 theDensityIdx = bld->GetCoupleIndexes();
131 131
132 scTracks.reserve(10); 132 scTracks.reserve(10);
133 secParticles.reserve(12); 133 secParticles.reserve(12);
134 emModels = new std::vector<G4VEmModel*>; 134 emModels = new std::vector<G4VEmModel*>;
135 } 135 }
136 136
137 //....oooOO0OOooo........oooOO0OOooo........oo 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 138
139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 139 G4VEnergyLossProcess::~G4VEnergyLossProcess()
140 { 140 {
141 if (isMaster) { 141 if (isMaster) {
142 if(nullptr == baseParticle) { delete theDa 142 if(nullptr == baseParticle) { delete theData; }
143 delete theEnergyOfCrossSectionMax; 143 delete theEnergyOfCrossSectionMax;
144 if(nullptr != fXSpeaks) { 144 if(nullptr != fXSpeaks) {
145 for(auto const & v : *fXSpeaks) { delete 145 for(auto const & v : *fXSpeaks) { delete v; }
146 delete fXSpeaks; 146 delete fXSpeaks;
147 } 147 }
148 } 148 }
149 delete modelManager; 149 delete modelManager;
150 delete biasManager; 150 delete biasManager;
151 delete scoffRegions; 151 delete scoffRegions;
152 delete emModels; 152 delete emModels;
153 lManager->DeRegister(this); 153 lManager->DeRegister(this);
154 } 154 }
155 155
156 //....oooOO0OOooo........oooOO0OOooo........oo 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 157
158 G4double G4VEnergyLossProcess::MinPrimaryEnerg 158 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
159 159 const G4Material*,
160 160 G4double cut)
161 { 161 {
162 return cut; 162 return cut;
163 } 163 }
164 164
165 //....oooOO0OOooo........oooOO0OOooo........oo 165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 166
167 void G4VEnergyLossProcess::AddEmModel(G4int or 167 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* ptr,
168 G4VEmFlu 168 G4VEmFluctuationModel* fluc,
169 const G4 169 const G4Region* region)
170 { 170 {
171 if(nullptr == ptr) { return; } 171 if(nullptr == ptr) { return; }
172 G4VEmFluctuationModel* afluc = (nullptr == f 172 G4VEmFluctuationModel* afluc = (nullptr == fluc) ? fluctModel : fluc;
173 modelManager->AddEmModel(order, ptr, afluc, 173 modelManager->AddEmModel(order, ptr, afluc, region);
174 ptr->SetParticleChange(pParticleChange, aflu 174 ptr->SetParticleChange(pParticleChange, afluc);
175 } 175 }
176 176
177 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 178
179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod 179 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* ptr, G4int)
180 { 180 {
181 if(nullptr == ptr) { return; } 181 if(nullptr == ptr) { return; }
182 if(!emModels->empty()) { 182 if(!emModels->empty()) {
183 for(auto & em : *emModels) { if(em == ptr) 183 for(auto & em : *emModels) { if(em == ptr) { return; } }
184 } 184 }
185 emModels->push_back(ptr); 185 emModels->push_back(ptr);
186 } 186 }
187 187
188 //....oooOO0OOooo........oooOO0OOooo........oo 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 189
190 void G4VEnergyLossProcess::SetDynamicMassCharg 190 void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
191 191 G4double charge2ratio)
192 { 192 {
193 massRatio = massratio; 193 massRatio = massratio;
194 logMassRatio = G4Log(massRatio); 194 logMassRatio = G4Log(massRatio);
195 fFactor = charge2ratio*biasFactor; 195 fFactor = charge2ratio*biasFactor;
196 if(baseMat) { fFactor *= (*theDensityFactor) 196 if(baseMat) { fFactor *= (*theDensityFactor)[currentCoupleIndex]; }
197 chargeSqRatio = charge2ratio; 197 chargeSqRatio = charge2ratio;
198 reduceFactor = 1.0/(fFactor*massRatio); 198 reduceFactor = 1.0/(fFactor*massRatio);
199 } 199 }
200 200
201 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 202
203 void 203 void
204 G4VEnergyLossProcess::PreparePhysicsTable(cons 204 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
205 { 205 {
206 particle = G4EmTableUtil::CheckIon(this, &pa 206 particle = G4EmTableUtil::CheckIon(this, &part, particle,
207 verboseLe 207 verboseLevel, isIon);
208 208
209 if( particle != &part ) { 209 if( particle != &part ) {
210 if(!isIon) { lManager->RegisterExtraPartic 210 if(!isIon) { lManager->RegisterExtraParticle(&part, this); }
211 if(1 < verboseLevel) { 211 if(1 < verboseLevel) {
212 G4cout << "### G4VEnergyLossProcess::Pre 212 G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()"
213 << " interrupted for " << GetProc << 213 << " interrupted for "
214 << part.GetParticleName() << " is << 214 << part.GetParticleName() << " isIon=" << isIon << G4endl;
215 << " spline=" << spline << G4endl <<
216 } 215 }
217 return; 216 return;
218 } 217 }
219 218
220 tablesAreBuilt = false; 219 tablesAreBuilt = false;
221 if (GetProcessSubType() == fIonisation) { Se <<
222 220
223 G4LossTableBuilder* bld = lManager->GetTable 221 G4LossTableBuilder* bld = lManager->GetTableBuilder();
224 lManager->PreparePhysicsTable(&part, this); << 222 lManager->PreparePhysicsTable(&part, this, isMaster);
225 223
226 // Base particle and set of models can be de 224 // Base particle and set of models can be defined here
227 InitialiseEnergyLossProcess(particle, basePa 225 InitialiseEnergyLossProcess(particle, baseParticle);
228 226
229 // parameters of the process 227 // parameters of the process
230 if(!actLossFluc) { lossFluctuationFlag = the 228 if(!actLossFluc) { lossFluctuationFlag = theParameters->LossFluctuation(); }
231 useCutAsFinalRange = theParameters->UseCutAs << 229 rndmStepFlag = theParameters->UseCutAsFinalRange();
232 if(!actMinKinEnergy) { minKinEnergy = thePar 230 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
233 if(!actMaxKinEnergy) { maxKinEnergy = thePar 231 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
234 if(!actBinning) { nBins = theParameters->Num 232 if(!actBinning) { nBins = theParameters->NumberOfBins(); }
235 maxKinEnergyCSDA = theParameters->MaxEnergyF 233 maxKinEnergyCSDA = theParameters->MaxEnergyForCSDARange();
236 nBinsCSDA = theParameters->NumberOfBinsPerDe 234 nBinsCSDA = theParameters->NumberOfBinsPerDecade()
237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi 235 *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy));
238 if(!actLinLossLimit) { linLossLimit = thePar 236 if(!actLinLossLimit) { linLossLimit = theParameters->LinearLossLimit(); }
239 lambdaFactor = theParameters->LambdaFactor() 237 lambdaFactor = theParameters->LambdaFactor();
240 invLambdaFactor = 1.0/lambdaFactor; 238 invLambdaFactor = 1.0/lambdaFactor;
241 if(isMaster) { SetVerboseLevel(theParameters 239 if(isMaster) { SetVerboseLevel(theParameters->Verbose()); }
242 else { SetVerboseLevel(theParameters->Worker 240 else { SetVerboseLevel(theParameters->WorkerVerbose()); }
243 // integral option may be disabled 241 // integral option may be disabled
244 if(!theParameters->Integral()) { fXSType = f 242 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; }
245 243
246 theParameters->DefineRegParamForLoss(this); 244 theParameters->DefineRegParamForLoss(this);
247 245
248 fRangeEnergy = 0.0; 246 fRangeEnergy = 0.0;
249 247
250 G4double initialCharge = particle->GetPDGCha 248 G4double initialCharge = particle->GetPDGCharge();
251 G4double initialMass = particle->GetPDGMas 249 G4double initialMass = particle->GetPDGMass();
252 250
253 theParameters->FillStepFunction(particle, th 251 theParameters->FillStepFunction(particle, this);
254 252
255 // parameters for scaling from the base part 253 // parameters for scaling from the base particle
256 if (nullptr != baseParticle) { 254 if (nullptr != baseParticle) {
257 massRatio = (baseParticle->GetPDGMass() 255 massRatio = (baseParticle->GetPDGMass())/initialMass;
258 logMassRatio = G4Log(massRatio); 256 logMassRatio = G4Log(massRatio);
259 G4double q = initialCharge/baseParticle->G 257 G4double q = initialCharge/baseParticle->GetPDGCharge();
260 chargeSqRatio = q*q; 258 chargeSqRatio = q*q;
261 if(chargeSqRatio > 0.0) { reduceFactor = 1 259 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
262 } 260 }
263 lowestKinEnergy = (initialMass < CLHEP::MeV) 261 lowestKinEnergy = (initialMass < CLHEP::MeV)
264 ? theParameters->LowestElectronEnergy() 262 ? theParameters->LowestElectronEnergy()
265 : theParameters->LowestMuHadEnergy(); 263 : theParameters->LowestMuHadEnergy();
266 264
267 // Tables preparation 265 // Tables preparation
268 if (isMaster && nullptr == baseParticle) { 266 if (isMaster && nullptr == baseParticle) {
269 if(nullptr == theData) { theData = new G4E 267 if(nullptr == theData) { theData = new G4EmDataHandler(7); }
270 268
271 if(nullptr != theDEDXTable && isIonisation 269 if(nullptr != theDEDXTable && isIonisation) {
272 if(nullptr != theIonisationTable && theD 270 if(nullptr != theIonisationTable && theDEDXTable != theIonisationTable) {
273 theData->CleanTable(0); 271 theData->CleanTable(0);
274 theDEDXTable = theIonisationTable; 272 theDEDXTable = theIonisationTable;
275 theIonisationTable = nullptr; 273 theIonisationTable = nullptr;
276 } 274 }
277 } 275 }
278 276
279 theDEDXTable = theData->MakeTable(theDEDXT 277 theDEDXTable = theData->MakeTable(theDEDXTable, 0);
280 bld->InitialiseBaseMaterials(theDEDXTable) 278 bld->InitialiseBaseMaterials(theDEDXTable);
281 theData->UpdateTable(theIonisationTable, 1 279 theData->UpdateTable(theIonisationTable, 1);
282 280
283 if (theParameters->BuildCSDARange()) { 281 if (theParameters->BuildCSDARange()) {
284 theDEDXunRestrictedTable = theData->Make 282 theDEDXunRestrictedTable = theData->MakeTable(2);
285 if(isIonisation) { theCSDARangeTable = t 283 if(isIonisation) { theCSDARangeTable = theData->MakeTable(3); }
286 } 284 }
287 285
288 theLambdaTable = theData->MakeTable(4); 286 theLambdaTable = theData->MakeTable(4);
289 if(isIonisation) { 287 if(isIonisation) {
290 theRangeTableForLoss = theData->MakeTabl 288 theRangeTableForLoss = theData->MakeTable(5);
291 theInverseRangeTable = theData->MakeTabl 289 theInverseRangeTable = theData->MakeTable(6);
292 } 290 }
293 } 291 }
294 292
295 // forced biasing 293 // forced biasing
296 if(nullptr != biasManager) { 294 if(nullptr != biasManager) {
297 biasManager->Initialise(part,GetProcessNam 295 biasManager->Initialise(part,GetProcessName(),verboseLevel);
298 biasFlag = false; 296 biasFlag = false;
299 } 297 }
300 baseMat = bld->GetBaseMaterialFlag(); 298 baseMat = bld->GetBaseMaterialFlag();
301 numberOfModels = modelManager->NumberOfModel 299 numberOfModels = modelManager->NumberOfModels();
302 currentModel = modelManager->GetModel(0); 300 currentModel = modelManager->GetModel(0);
303 G4EmTableUtil::UpdateModels(this, modelManag 301 G4EmTableUtil::UpdateModels(this, modelManager, maxKinEnergy,
304 numberOfModels, 302 numberOfModels, secID, biasID,
305 mainSecondaries, 303 mainSecondaries, baseMat, isMaster,
306 theParameters->U 304 theParameters->UseAngularGeneratorForIonisation());
307 theCuts = modelManager->Initialise(particle, 305 theCuts = modelManager->Initialise(particle, secondaryParticle,
308 verboseLe 306 verboseLevel);
309 // subcut processor 307 // subcut processor
310 if(isIonisation) { 308 if(isIonisation) {
311 subcutProducer = lManager->SubCutProducer( 309 subcutProducer = lManager->SubCutProducer();
312 } 310 }
313 if(1 == nSCoffRegions) { 311 if(1 == nSCoffRegions) {
314 if((*scoffRegions)[0]->GetName() == "Defau 312 if((*scoffRegions)[0]->GetName() == "DefaultRegionForTheWorld") {
315 delete scoffRegions; 313 delete scoffRegions;
316 scoffRegions = nullptr; 314 scoffRegions = nullptr;
317 nSCoffRegions = 0; 315 nSCoffRegions = 0;
318 } 316 }
319 } 317 }
320 318
321 if(1 < verboseLevel) { 319 if(1 < verboseLevel) {
322 G4cout << "G4VEnergyLossProcess::PrepearPh 320 G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done "
323 << " for " << GetProcessName() << " << 321 << " for local " << particle->GetParticleName()
324 << " isIon= " << isIon << " spline= << 322 << " isIon= " << isIon;
325 if(baseParticle) { 323 if(baseParticle) {
326 G4cout << "; base: " << baseParticle->Ge 324 G4cout << "; base: " << baseParticle->GetParticleName();
327 } 325 }
328 G4cout << G4endl; <<
329 G4cout << " chargeSqRatio= " << chargeSqRa 326 G4cout << " chargeSqRatio= " << chargeSqRatio
330 << " massRatio= " << massRatio 327 << " massRatio= " << massRatio
331 << " reduceFactor= " << reduceFacto 328 << " reduceFactor= " << reduceFactor << G4endl;
332 if (nSCoffRegions > 0) { 329 if (nSCoffRegions > 0) {
333 G4cout << " SubCut secondary production 330 G4cout << " SubCut secondary production is ON for regions: " << G4endl;
334 for (G4int i=0; i<nSCoffRegions; ++i) { 331 for (G4int i=0; i<nSCoffRegions; ++i) {
335 const G4Region* r = (*scoffRegions)[i] 332 const G4Region* r = (*scoffRegions)[i];
336 G4cout << " " << r->GetName( 333 G4cout << " " << r->GetName() << G4endl;
337 } 334 }
338 } else if(nullptr != subcutProducer) { 335 } else if(nullptr != subcutProducer) {
339 G4cout << " SubCut secondary production 336 G4cout << " SubCut secondary production is ON for all regions" << G4endl;
340 } 337 }
341 } 338 }
342 } 339 }
343 340
344 //....oooOO0OOooo........oooOO0OOooo........oo 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 342
346 void G4VEnergyLossProcess::BuildPhysicsTable(c 343 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
347 { 344 {
348 if(1 < verboseLevel) { 345 if(1 < verboseLevel) {
349 G4cout << "### G4VEnergyLossProcess::Build 346 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
350 << GetProcessName() 347 << GetProcessName()
351 << " and particle " << part.GetPart 348 << " and particle " << part.GetParticleName()
352 << "; the first particle " << parti << 349 << "; local: " << particle->GetParticleName();
353 if(baseParticle) { 350 if(baseParticle) {
354 G4cout << "; base: " << baseParticle->Ge 351 G4cout << "; base: " << baseParticle->GetParticleName();
355 } 352 }
356 G4cout << G4endl; << 353 G4cout << " TablesAreBuilt= " << tablesAreBuilt
357 G4cout << " TablesAreBuilt= " << tables << 354 << " isIon= " << isIon << " " << this << G4endl;
358 << " spline=" << spline << " ptr: " <<
359 } 355 }
360 356
361 if(&part == particle) { 357 if(&part == particle) {
>> 358
362 if(isMaster) { 359 if(isMaster) {
363 lManager->BuildPhysicsTable(particle, th 360 lManager->BuildPhysicsTable(particle, this);
364 361
365 } else { 362 } else {
366 const auto masterProcess = 363 const auto masterProcess =
367 static_cast<const G4VEnergyLossProcess 364 static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
368 365
369 numberOfModels = modelManager->NumberOfM 366 numberOfModels = modelManager->NumberOfModels();
370 G4EmTableUtil::BuildLocalElossProcess(th 367 G4EmTableUtil::BuildLocalElossProcess(this, masterProcess,
371 pa 368 particle, numberOfModels);
372 tablesAreBuilt = true; 369 tablesAreBuilt = true;
373 baseMat = masterProcess->UseBaseMaterial 370 baseMat = masterProcess->UseBaseMaterial();
374 lManager->LocalPhysicsTables(particle, t 371 lManager->LocalPhysicsTables(particle, this);
375 } 372 }
376 373
377 // needs to be done only once 374 // needs to be done only once
378 safetyHelper->InitialiseHelper(); 375 safetyHelper->InitialiseHelper();
379 } 376 }
380 // Added tracking cut to avoid tracking arti 377 // Added tracking cut to avoid tracking artifacts
381 // and identified deexcitation flag 378 // and identified deexcitation flag
382 if(isIonisation) { 379 if(isIonisation) {
383 atomDeexcitation = lManager->AtomDeexcitat 380 atomDeexcitation = lManager->AtomDeexcitation();
384 if(nullptr != atomDeexcitation) { 381 if(nullptr != atomDeexcitation) {
385 if(atomDeexcitation->IsPIXEActive()) { u 382 if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; }
386 } 383 }
387 } 384 }
388 385
389 // protection against double printout 386 // protection against double printout
390 if(theParameters->IsPrintLocked()) { return; 387 if(theParameters->IsPrintLocked()) { return; }
391 388
392 // explicitly defined printout by particle n 389 // explicitly defined printout by particle name
393 G4String num = part.GetParticleName(); 390 G4String num = part.GetParticleName();
394 if(1 < verboseLevel || 391 if(1 < verboseLevel ||
395 (0 < verboseLevel && (num == "e-" || 392 (0 < verboseLevel && (num == "e-" ||
396 num == "e+" || n 393 num == "e+" || num == "mu+" ||
397 num == "mu-" || n 394 num == "mu-" || num == "proton"||
398 num == "pi+" || n 395 num == "pi+" || num == "pi-" ||
399 num == "kaon+" || n 396 num == "kaon+" || num == "kaon-" ||
400 num == "alpha" || n 397 num == "alpha" || num == "anti_proton" ||
401 num == "GenericIon" 398 num == "GenericIon"|| num == "alpha+" ))) {
402 StreamInfo(G4cout, part); 399 StreamInfo(G4cout, part);
403 } 400 }
404 if(1 < verboseLevel) { 401 if(1 < verboseLevel) {
405 G4cout << "### G4VEnergyLossProcess::Build 402 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
406 << GetProcessName() 403 << GetProcessName()
407 << " and particle " << part.GetPart 404 << " and particle " << part.GetParticleName();
408 if(isIonisation) { G4cout << " isIonisati 405 if(isIonisation) { G4cout << " isIonisation flag=1"; }
409 G4cout << " baseMat=" << baseMat << G4endl 406 G4cout << " baseMat=" << baseMat << G4endl;
410 } 407 }
411 } 408 }
412 409
413 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414 411
415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED 412 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
416 { 413 {
417 G4PhysicsTable* table = nullptr; 414 G4PhysicsTable* table = nullptr;
418 G4double emax = maxKinEnergy; 415 G4double emax = maxKinEnergy;
419 G4int bin = nBins; 416 G4int bin = nBins;
420 417
421 if(fTotal == tType) { 418 if(fTotal == tType) {
422 emax = maxKinEnergyCSDA; 419 emax = maxKinEnergyCSDA;
423 bin = nBinsCSDA; 420 bin = nBinsCSDA;
424 table = theDEDXunRestrictedTable; 421 table = theDEDXunRestrictedTable;
425 } else if(fRestricted == tType) { 422 } else if(fRestricted == tType) {
426 table = theDEDXTable; 423 table = theDEDXTable;
427 } else { 424 } else {
428 G4cout << "G4VEnergyLossProcess::BuildDEDX 425 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
429 << tType << G4endl; 426 << tType << G4endl;
430 } 427 }
431 if(1 < verboseLevel) { 428 if(1 < verboseLevel) {
432 G4cout << "G4VEnergyLossProcess::BuildDEDX 429 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
433 << " for " << GetProcessName() 430 << " for " << GetProcessName()
434 << " and " << particle->GetParticle << 431 << " and " << particle->GetParticleName() << G4endl;
435 << "spline=" << spline << G4endl; <<
436 } 432 }
437 if(nullptr == table) { return table; } 433 if(nullptr == table) { return table; }
438 434
439 G4LossTableBuilder* bld = lManager->GetTable 435 G4LossTableBuilder* bld = lManager->GetTableBuilder();
440 G4EmTableUtil::BuildDEDXTable(this, particle 436 G4EmTableUtil::BuildDEDXTable(this, particle, modelManager, bld,
441 table, minKinE 437 table, minKinEnergy, emax, bin,
442 verboseLevel, 438 verboseLevel, tType, spline);
443 return table; 439 return table;
444 } 440 }
445 441
446 //....oooOO0OOooo........oooOO0OOooo........oo 442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 443
448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam 444 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType)
449 { 445 {
450 if(nullptr == theLambdaTable) { return theLa 446 if(nullptr == theLambdaTable) { return theLambdaTable; }
451 447
452 G4double scale = theParameters->MaxKinEnergy 448 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
453 G4int nbin = 449 G4int nbin =
454 theParameters->NumberOfBinsPerDecade()*G4l 450 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale));
455 scale = nbin/G4Log(scale); 451 scale = nbin/G4Log(scale);
456 452
457 G4LossTableBuilder* bld = lManager->GetTable 453 G4LossTableBuilder* bld = lManager->GetTableBuilder();
458 G4EmTableUtil::BuildLambdaTable(this, partic 454 G4EmTableUtil::BuildLambdaTable(this, particle, modelManager,
459 bld, theLamb 455 bld, theLambdaTable, theCuts,
460 minKinEnergy 456 minKinEnergy, maxKinEnergy, scale,
461 verboseLevel 457 verboseLevel, spline);
462 return theLambdaTable; 458 return theLambdaTable;
463 } 459 }
464 460
465 //....oooOO0OOooo........oooOO0OOooo........oo 461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 462
467 void G4VEnergyLossProcess::StreamInfo(std::ost 463 void G4VEnergyLossProcess::StreamInfo(std::ostream& out,
468 const G4ParticleDefinition& pa 464 const G4ParticleDefinition& part, G4bool rst) const
469 { 465 {
470 G4String indent = (rst ? " " : ""); 466 G4String indent = (rst ? " " : "");
471 out << std::setprecision(6); 467 out << std::setprecision(6);
472 out << G4endl << indent << GetProcessName() 468 out << G4endl << indent << GetProcessName() << ": ";
473 if (!rst) out << " for " << part.GetParticle 469 if (!rst) out << " for " << part.GetParticleName();
474 out << " XStype:" << fXSType 470 out << " XStype:" << fXSType
475 << " SubType=" << GetProcessSubType() < 471 << " SubType=" << GetProcessSubType() << G4endl
476 << " dE/dx and range tables from " 472 << " dE/dx and range tables from "
477 << G4BestUnit(minKinEnergy,"Energy") 473 << G4BestUnit(minKinEnergy,"Energy")
478 << " to " << G4BestUnit(maxKinEnergy,"En 474 << " to " << G4BestUnit(maxKinEnergy,"Energy")
479 << " in " << nBins << " bins" << G4endl 475 << " in " << nBins << " bins" << G4endl
480 << " Lambda tables from threshold t 476 << " Lambda tables from threshold to "
481 << G4BestUnit(maxKinEnergy,"Energy") 477 << G4BestUnit(maxKinEnergy,"Energy")
482 << ", " << theParameters->NumberOfBinsPe 478 << ", " << theParameters->NumberOfBinsPerDecade()
483 << " bins/decade, spline: " << spline 479 << " bins/decade, spline: " << spline
484 << G4endl; 480 << G4endl;
485 if(nullptr != theRangeTableForLoss && isIoni 481 if(nullptr != theRangeTableForLoss && isIonisation) {
486 out << " StepFunction=(" << dRoverRan 482 out << " StepFunction=(" << dRoverRange << ", "
487 << finalRange/mm << " mm)" 483 << finalRange/mm << " mm)"
488 << ", integ: " << fXSType 484 << ", integ: " << fXSType
489 << ", fluct: " << lossFluctuationFlag 485 << ", fluct: " << lossFluctuationFlag
490 << ", linLossLim= " << linLossLimit 486 << ", linLossLim= " << linLossLimit
491 << G4endl; 487 << G4endl;
492 } 488 }
493 StreamProcessInfo(out); 489 StreamProcessInfo(out);
494 modelManager->DumpModelList(out, verboseLeve 490 modelManager->DumpModelList(out, verboseLevel);
495 if(nullptr != theCSDARangeTable && isIonisat 491 if(nullptr != theCSDARangeTable && isIonisation) {
496 out << " CSDA range table up" 492 out << " CSDA range table up"
497 << " to " << G4BestUnit(maxKinEnergyCS 493 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
498 << " in " << nBinsCSDA << " bins" << G 494 << " in " << nBinsCSDA << " bins" << G4endl;
499 } 495 }
500 if(nSCoffRegions>0 && isIonisation) { 496 if(nSCoffRegions>0 && isIonisation) {
501 out << " Subcutoff sampling in " << n 497 out << " Subcutoff sampling in " << nSCoffRegions
502 << " regions" << G4endl; 498 << " regions" << G4endl;
503 } 499 }
504 if(2 < verboseLevel) { 500 if(2 < verboseLevel) {
505 for(std::size_t i=0; i<7; ++i) { 501 for(std::size_t i=0; i<7; ++i) {
506 auto ta = theData->Table(i); 502 auto ta = theData->Table(i);
507 out << " " << tnames[i] << " addres 503 out << " " << tnames[i] << " address: " << ta << G4endl;
508 if(nullptr != ta) { out << *ta << G4endl 504 if(nullptr != ta) { out << *ta << G4endl; }
509 } 505 }
510 } 506 }
511 } 507 }
512 508
513 //....oooOO0OOooo........oooOO0OOooo........oo 509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
514 510
515 void G4VEnergyLossProcess::ActivateSubCutoff(c 511 void G4VEnergyLossProcess::ActivateSubCutoff(const G4Region* r)
516 { 512 {
517 if(nullptr == scoffRegions) { 513 if(nullptr == scoffRegions) {
518 scoffRegions = new std::vector<const G4Reg 514 scoffRegions = new std::vector<const G4Region*>;
519 } 515 }
520 // the region is in the list 516 // the region is in the list
521 if(!scoffRegions->empty()) { 517 if(!scoffRegions->empty()) {
522 for (auto & reg : *scoffRegions) { 518 for (auto & reg : *scoffRegions) {
523 if (reg == r) { return; } 519 if (reg == r) { return; }
524 } 520 }
525 } 521 }
526 // new region 522 // new region
527 scoffRegions->push_back(r); 523 scoffRegions->push_back(r);
528 ++nSCoffRegions; 524 ++nSCoffRegions;
529 } 525 }
530 526
531 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
532 528
533 G4bool G4VEnergyLossProcess::IsRegionForCubcut 529 G4bool G4VEnergyLossProcess::IsRegionForCubcutProcessor(const G4Track& aTrack)
534 { 530 {
535 if(0 == nSCoffRegions) { return true; } 531 if(0 == nSCoffRegions) { return true; }
536 const G4Region* r = aTrack.GetVolume()->GetL 532 const G4Region* r = aTrack.GetVolume()->GetLogicalVolume()->GetRegion();
537 for(auto & reg : *scoffRegions) { 533 for(auto & reg : *scoffRegions) {
538 if(r == reg) { return true; } 534 if(r == reg) { return true; }
539 } 535 }
540 return false; 536 return false;
541 } 537 }
542 538
543 //....oooOO0OOooo........oooOO0OOooo........oo 539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544 540
545 void G4VEnergyLossProcess::StartTracking(G4Tra 541 void G4VEnergyLossProcess::StartTracking(G4Track* track)
546 { 542 {
547 // reset parameters for the new track 543 // reset parameters for the new track
548 theNumberOfInteractionLengthLeft = -1.0; 544 theNumberOfInteractionLengthLeft = -1.0;
549 mfpKinEnergy = DBL_MAX; 545 mfpKinEnergy = DBL_MAX;
550 preStepLambda = 0.0; <<
551 currentCouple = nullptr; 546 currentCouple = nullptr;
552 547
553 // reset ion 548 // reset ion
554 if(isIon) { 549 if(isIon) {
555 const G4double newmass = track->GetDefinit 550 const G4double newmass = track->GetDefinition()->GetPDGMass();
556 massRatio = (nullptr == baseParticle) ? CL << 551 if(nullptr != baseParticle) {
557 : baseParticle->GetPDGMass()/newmass; << 552 massRatio = baseParticle->GetPDGMass()/newmass;
558 logMassRatio = G4Log(massRatio); << 553 logMassRatio = G4Log(massRatio);
>> 554 } else if(isIon) {
>> 555 massRatio = CLHEP::proton_mass_c2/newmass;
>> 556 logMassRatio = G4Log(massRatio);
>> 557 } else {
>> 558 massRatio = 1.0;
>> 559 logMassRatio = 0.0;
>> 560 }
559 } 561 }
560 // forced biasing only for primary particles 562 // forced biasing only for primary particles
561 if(nullptr != biasManager) { 563 if(nullptr != biasManager) {
562 if(0 == track->GetParentID()) { 564 if(0 == track->GetParentID()) {
563 biasFlag = true; 565 biasFlag = true;
564 biasManager->ResetForcedInteraction(); 566 biasManager->ResetForcedInteraction();
565 } 567 }
566 } 568 }
567 } 569 }
568 570
569 //....oooOO0OOooo........oooOO0OOooo........oo 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
570 572
571 G4double G4VEnergyLossProcess::AlongStepGetPhy 573 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
572 const G4Track& tr << 574 const G4Track&,G4double,G4double,G4double&,
573 G4GPILSelection* 575 G4GPILSelection* selection)
574 { 576 {
575 G4double x = DBL_MAX; 577 G4double x = DBL_MAX;
576 *selection = aGPILSelection; 578 *selection = aGPILSelection;
577 if(isIonisation && currentModel->IsActive(pr 579 if(isIonisation && currentModel->IsActive(preStepScaledEnergy)) {
578 GetScaledRangeForScaledEnergy(preStepScale << 580 GetScaledRangeForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy);
579 x = (useCutAsFinalRange) ? std::min(finalR << 581 const G4double finR = (rndmStepFlag) ? std::min(finalRange,
580 currentCouple->GetProductionCuts()->GetP 582 currentCouple->GetProductionCuts()->GetProductionCut(1)) : finalRange;
581 x = (fRange > x) ? fRange*dRoverRange + x* << 583 x = (fRange > finR) ?
582 : fRange; << 584 fRange*dRoverRange + finR*(1.0-dRoverRange)*(2.0-finR/fRange) : fRange;
583 /* <<
584 G4cout<<"AlongStepGPIL: " << GetProcessN <<
585 << " fRange=" << fRange << " finR=" << finR <<
586 */ <<
587 } 585 }
>> 586 //G4cout<<"AlongStepGPIL: " << GetProcessName()<<": e= "<<preStepKinEnergy
>> 587 //<<" stepLimit= "<<x<<G4endl;
588 return x; 588 return x;
589 } 589 }
590 590
591 //....oooOO0OOooo........oooOO0OOooo........oo 591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
592 592
593 G4double G4VEnergyLossProcess::PostStepGetPhys 593 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
594 const G4Track& tr 594 const G4Track& track,
595 G4double previo 595 G4double previousStepSize,
596 G4ForceCondition* 596 G4ForceCondition* condition)
597 { 597 {
598 // condition is set to "Not Forced" 598 // condition is set to "Not Forced"
599 *condition = NotForced; 599 *condition = NotForced;
600 G4double x = DBL_MAX; 600 G4double x = DBL_MAX;
601 601
602 // initialisation of material, mass, charge, 602 // initialisation of material, mass, charge, model
603 // at the beginning of the step 603 // at the beginning of the step
604 DefineMaterial(track.GetMaterialCutsCouple() 604 DefineMaterial(track.GetMaterialCutsCouple());
605 preStepKinEnergy = track.GetKineticEne 605 preStepKinEnergy = track.GetKineticEnergy();
>> 606 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
606 preStepScaledEnergy = preStepKinEnergy*ma 607 preStepScaledEnergy = preStepKinEnergy*massRatio;
>> 608 preStepLogScaledEnergy = preStepLogKinEnergy + logMassRatio;
607 SelectModel(preStepScaledEnergy); 609 SelectModel(preStepScaledEnergy);
608 610
609 if(!currentModel->IsActive(preStepScaledEner 611 if(!currentModel->IsActive(preStepScaledEnergy)) {
610 theNumberOfInteractionLengthLeft = -1.0; 612 theNumberOfInteractionLengthLeft = -1.0;
611 mfpKinEnergy = DBL_MAX; <<
612 preStepLambda = 0.0; <<
613 currentInteractionLength = DBL_MAX; 613 currentInteractionLength = DBL_MAX;
614 return x; << 614 return x;
615 } 615 }
616 616
617 // change effective charge of a charged part 617 // change effective charge of a charged particle on fly
618 if(isIon) { 618 if(isIon) {
619 const G4double q2 = currentModel->ChargeSq 619 const G4double q2 = currentModel->ChargeSquareRatio(track);
620 fFactor = q2*biasFactor; << 620 if(q2 != chargeSqRatio) {
621 if(baseMat) { fFactor *= (*theDensityFacto << 621 fFactor *= q2/chargeSqRatio;
622 reduceFactor = 1.0/(fFactor*massRatio); << 622 reduceFactor = 1.0/(fFactor*massRatio);
>> 623 chargeSqRatio = q2;
>> 624 }
623 if (lossFluctuationFlag) { 625 if (lossFluctuationFlag) {
624 auto fluc = currentModel->GetModelOfFluc 626 auto fluc = currentModel->GetModelOfFluctuations();
625 fluc->SetParticleAndCharge(track.GetDefi 627 fluc->SetParticleAndCharge(track.GetDefinition(), q2);
626 } 628 }
627 } 629 }
628 630
629 // forced biasing only for primary particles 631 // forced biasing only for primary particles
630 if(biasManager) { 632 if(biasManager) {
631 if(0 == track.GetParentID() && biasFlag && 633 if(0 == track.GetParentID() && biasFlag &&
632 biasManager->ForcedInteractionRegion((G 634 biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
633 return biasManager->GetStepLimit((G4int) 635 return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize);
634 } 636 }
635 } 637 }
636 638
637 ComputeLambdaForScaledEnergy(preStepScaledEn << 639 // compute mean free path
638 << 640 ComputeLambdaForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy);
>> 641
639 // zero cross section 642 // zero cross section
640 if(preStepLambda <= 0.0) { 643 if(preStepLambda <= 0.0) {
641 theNumberOfInteractionLengthLeft = -1.0; 644 theNumberOfInteractionLengthLeft = -1.0;
642 currentInteractionLength = DBL_MAX; 645 currentInteractionLength = DBL_MAX;
643 } else { 646 } else {
644 647
645 // non-zero cross section 648 // non-zero cross section
646 if (theNumberOfInteractionLengthLeft < 0.0 649 if (theNumberOfInteractionLengthLeft < 0.0) {
647 650
648 // beggining of tracking (or just after 651 // beggining of tracking (or just after DoIt of this process)
649 theNumberOfInteractionLengthLeft = -G4Lo 652 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
650 theInitialNumberOfInteractionLength = th 653 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
651 654
652 } else if(currentInteractionLength < DBL_M 655 } else if(currentInteractionLength < DBL_MAX) {
653 656
654 // subtract NumberOfInteractionLengthLef 657 // subtract NumberOfInteractionLengthLeft using previous step
655 theNumberOfInteractionLengthLeft -= 658 theNumberOfInteractionLengthLeft -=
656 previousStepSize/currentInteractionLen 659 previousStepSize/currentInteractionLength;
657 660
658 theNumberOfInteractionLengthLeft = 661 theNumberOfInteractionLengthLeft =
659 std::max(theNumberOfInteractionLengthL 662 std::max(theNumberOfInteractionLengthLeft, 0.0);
660 } 663 }
661 664
662 // new mean free path and step limit 665 // new mean free path and step limit
663 currentInteractionLength = 1.0/preStepLamb 666 currentInteractionLength = 1.0/preStepLambda;
664 x = theNumberOfInteractionLengthLeft * cur 667 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
665 } 668 }
666 #ifdef G4VERBOSE 669 #ifdef G4VERBOSE
667 if (verboseLevel>2) { 670 if (verboseLevel>2) {
668 G4cout << "G4VEnergyLossProcess::PostStepG 671 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
669 G4cout << "[ " << GetProcessName() << "]" 672 G4cout << "[ " << GetProcessName() << "]" << G4endl;
670 G4cout << " for " << track.GetDefinition() 673 G4cout << " for " << track.GetDefinition()->GetParticleName()
671 << " in Material " << currentMate 674 << " in Material " << currentMaterial->GetName()
672 << " Ekin(MeV)= " << preStepKinEner 675 << " Ekin(MeV)= " << preStepKinEnergy/MeV
673 << " track material: " << track.Get 676 << " track material: " << track.GetMaterial()->GetName()
674 <<G4endl; 677 <<G4endl;
675 G4cout << "MeanFreePath = " << currentInte 678 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
676 << "InteractionLength= " << x/cm << 679 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
677 } 680 }
678 #endif 681 #endif
679 return x; 682 return x;
680 } 683 }
681 684
682 //....oooOO0OOooo........oooOO0OOooo........oo 685 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683 686
684 void 687 void
685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 688 G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e, G4double loge)
686 { 689 {
687 // cross section increased with energy 690 // cross section increased with energy
688 if(fXSType == fEmIncreasing) { 691 if(fXSType == fEmIncreasing) {
689 if(e*invLambdaFactor < mfpKinEnergy) { 692 if(e*invLambdaFactor < mfpKinEnergy) {
690 preStepLambda = GetLambdaForScaledEnergy << 693 mfpKinEnergy = e;
691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 694 preStepLambda = GetLambdaForScaledEnergy(e, loge);
692 } 695 }
693 696
694 // cross section has one peak 697 // cross section has one peak
695 } else if(fXSType == fEmOnePeak) { 698 } else if(fXSType == fEmOnePeak) {
696 const G4double epeak = (*theEnergyOfCrossS 699 const G4double epeak = (*theEnergyOfCrossSectionMax)[basedCoupleIndex];
697 if(e <= epeak) { 700 if(e <= epeak) {
698 if(e*invLambdaFactor < mfpKinEnergy) { 701 if(e*invLambdaFactor < mfpKinEnergy) {
699 preStepLambda = GetLambdaForScaledEner << 702 mfpKinEnergy = e;
700 mfpKinEnergy = (preStepLambda > 0.0) ? << 703 preStepLambda = GetLambdaForScaledEnergy(e, loge);
701 } 704 }
702 } else if(e < mfpKinEnergy) { 705 } else if(e < mfpKinEnergy) {
703 const G4double e1 = std::max(epeak, e*la 706 const G4double e1 = std::max(epeak, e*lambdaFactor);
704 mfpKinEnergy = e1; 707 mfpKinEnergy = e1;
705 preStepLambda = GetLambdaForScaledEnergy 708 preStepLambda = GetLambdaForScaledEnergy(e1);
706 } 709 }
707 710
708 // cross section has more than one peaks 711 // cross section has more than one peaks
709 } else if(fXSType == fEmTwoPeaks) { 712 } else if(fXSType == fEmTwoPeaks) {
710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple 713 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCoupleIndex];
711 const G4double e1peak = xs->e1peak; 714 const G4double e1peak = xs->e1peak;
712 715
713 // below the 1st peak 716 // below the 1st peak
714 if(e <= e1peak) { 717 if(e <= e1peak) {
715 if(e*invLambdaFactor < mfpKinEnergy) { 718 if(e*invLambdaFactor < mfpKinEnergy) {
716 preStepLambda = GetLambdaForScaledEner << 719 mfpKinEnergy = e;
717 mfpKinEnergy = (preStepLambda > 0.0) ? << 720 preStepLambda = GetLambdaForScaledEnergy(e, loge);
718 } 721 }
719 return; 722 return;
720 } 723 }
721 const G4double e1deep = xs->e1deep; 724 const G4double e1deep = xs->e1deep;
722 // above the 1st peak, below the deep 725 // above the 1st peak, below the deep
723 if(e <= e1deep) { 726 if(e <= e1deep) {
724 if(mfpKinEnergy >= e1deep || e <= mfpKin 727 if(mfpKinEnergy >= e1deep || e <= mfpKinEnergy) {
725 const G4double e1 = std::max(e1peak, e 728 const G4double e1 = std::max(e1peak, e*lambdaFactor);
726 mfpKinEnergy = e1; 729 mfpKinEnergy = e1;
727 preStepLambda = GetLambdaForScaledEner 730 preStepLambda = GetLambdaForScaledEnergy(e1);
728 } 731 }
729 return; 732 return;
730 } 733 }
731 const G4double e2peak = xs->e2peak; 734 const G4double e2peak = xs->e2peak;
732 // above the deep, below 2nd peak 735 // above the deep, below 2nd peak
733 if(e <= e2peak) { 736 if(e <= e2peak) {
734 if(e*invLambdaFactor < mfpKinEnergy) { 737 if(e*invLambdaFactor < mfpKinEnergy) {
735 mfpKinEnergy = e; 738 mfpKinEnergy = e;
736 preStepLambda = GetLambdaForScaledEner << 739 preStepLambda = GetLambdaForScaledEnergy(e, loge);
737 } 740 }
738 return; 741 return;
739 } 742 }
740 const G4double e2deep = xs->e2deep; 743 const G4double e2deep = xs->e2deep;
741 // above the 2nd peak, below the deep 744 // above the 2nd peak, below the deep
742 if(e <= e2deep) { 745 if(e <= e2deep) {
743 if(mfpKinEnergy >= e2deep || e <= mfpKin 746 if(mfpKinEnergy >= e2deep || e <= mfpKinEnergy) {
744 const G4double e1 = std::max(e2peak, e 747 const G4double e1 = std::max(e2peak, e*lambdaFactor);
745 mfpKinEnergy = e1; 748 mfpKinEnergy = e1;
746 preStepLambda = GetLambdaForScaledEner 749 preStepLambda = GetLambdaForScaledEnergy(e1);
747 } 750 }
748 return; 751 return;
749 } 752 }
750 const G4double e3peak = xs->e3peak; 753 const G4double e3peak = xs->e3peak;
751 // above the deep, below 3d peak 754 // above the deep, below 3d peak
752 if(e <= e3peak) { 755 if(e <= e3peak) {
753 if(e*invLambdaFactor < mfpKinEnergy) { 756 if(e*invLambdaFactor < mfpKinEnergy) {
754 mfpKinEnergy = e; 757 mfpKinEnergy = e;
755 preStepLambda = GetLambdaForScaledEner << 758 preStepLambda = GetLambdaForScaledEnergy(e, loge);
756 } 759 }
757 return; 760 return;
758 } 761 }
759 // above 3d peak 762 // above 3d peak
760 if(e <= mfpKinEnergy) { 763 if(e <= mfpKinEnergy) {
761 const G4double e1 = std::max(e3peak, e*l 764 const G4double e1 = std::max(e3peak, e*lambdaFactor);
762 mfpKinEnergy = e1; 765 mfpKinEnergy = e1;
763 preStepLambda = GetLambdaForScaledEnergy 766 preStepLambda = GetLambdaForScaledEnergy(e1);
764 } 767 }
765 // integral method is not used 768 // integral method is not used
766 } else { 769 } else {
767 preStepLambda = GetLambdaForScaledEnergy(e << 770 preStepLambda = GetLambdaForScaledEnergy(e, loge);
768 } 771 }
769 } 772 }
770 773
771 //....oooOO0OOooo........oooOO0OOooo........oo 774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
772 775
773 G4VParticleChange* G4VEnergyLossProcess::Along 776 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
774 777 const G4Step& step)
775 { 778 {
776 fParticleChange.InitializeForAlongStep(track 779 fParticleChange.InitializeForAlongStep(track);
777 // The process has range table - calculate e 780 // The process has range table - calculate energy loss
778 if(!isIonisation || !currentModel->IsActive( 781 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) {
779 return &fParticleChange; 782 return &fParticleChange;
780 } 783 }
781 784
>> 785 // Get the actual (true) Step length
782 G4double length = step.GetStepLength(); 786 G4double length = step.GetStepLength();
>> 787 if(length <= 0.0) { return &fParticleChange; }
783 G4double eloss = 0.0; 788 G4double eloss = 0.0;
784 789
785 /* << 790 /*
786 if(-1 < verboseLevel) { 791 if(-1 < verboseLevel) {
787 const G4ParticleDefinition* d = track.GetP 792 const G4ParticleDefinition* d = track.GetParticleDefinition();
788 G4cout << "AlongStepDoIt for " 793 G4cout << "AlongStepDoIt for "
789 << GetProcessName() << " and partic 794 << GetProcessName() << " and particle " << d->GetParticleName()
790 << " eScaled(MeV)=" << preStepScal 795 << " eScaled(MeV)=" << preStepScaledEnergy/MeV
791 << " range(mm)=" << fRange/mm << " 796 << " range(mm)=" << fRange/mm << " s(mm)=" << length/mm
792 << " rf=" << reduceFactor << " q^ 797 << " rf=" << reduceFactor << " q^2=" << chargeSqRatio
793 << " md=" << d->GetPDGMass() << " 798 << " md=" << d->GetPDGMass() << " status=" << track.GetTrackStatus()
794 << " " << track.GetMaterial()->Get 799 << " " << track.GetMaterial()->GetName() << G4endl;
795 } 800 }
796 */ 801 */
797 const G4DynamicParticle* dynParticle = track 802 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
798 803
799 // define new weight for primary and seconda 804 // define new weight for primary and secondaries
800 G4double weight = fParticleChange.GetParentW 805 G4double weight = fParticleChange.GetParentWeight();
801 if(weightFlag) { 806 if(weightFlag) {
802 weight /= biasFactor; 807 weight /= biasFactor;
803 fParticleChange.ProposeWeight(weight); 808 fParticleChange.ProposeWeight(weight);
804 } 809 }
805 810
806 // stopping, check actual range and kinetic << 811 // stopping
807 if (length >= fRange || preStepKinEnergy <= 812 if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) {
808 eloss = preStepKinEnergy; 813 eloss = preStepKinEnergy;
809 if (useDeexcitation) { 814 if (useDeexcitation) {
810 atomDeexcitation->AlongStepDeexcitation( 815 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
811 816 eloss, (G4int)currentCoupleIndex);
812 if(scTracks.size() > 0) { FillSecondarie 817 if(scTracks.size() > 0) { FillSecondariesAlongStep(weight); }
813 eloss = std::max(eloss, 0.0); 818 eloss = std::max(eloss, 0.0);
814 } 819 }
815 fParticleChange.SetProposedKineticEnergy(0 820 fParticleChange.SetProposedKineticEnergy(0.0);
816 fParticleChange.ProposeLocalEnergyDeposit( 821 fParticleChange.ProposeLocalEnergyDeposit(eloss);
817 return &fParticleChange; 822 return &fParticleChange;
818 } 823 }
819 // zero step length with non-zero range <<
820 if(length <= 0.0) { return &fParticleChange; <<
821 <<
822 // Short step 824 // Short step
823 eloss = length*GetDEDXForScaledEnergy(preSte 825 eloss = length*GetDEDXForScaledEnergy(preStepScaledEnergy,
824 LogSca << 826 preStepLogScaledEnergy);
825 /* << 827 //G4cout << "Short STEP: eloss= " << eloss << G4endl;
826 G4cout << "##### Short STEP: eloss= " << elo << 828
827 << " Escaled=" << preStepScaledEnergy <<
828 << " R=" << fRange <<
829 << " L=" << length <<
830 << " fFactor=" << fFactor << " minE=" << mi <<
831 << " idxBase=" << basedCoupleIndex << G4end <<
832 */ <<
833 // Long step 829 // Long step
834 if(eloss > preStepKinEnergy*linLossLimit) { 830 if(eloss > preStepKinEnergy*linLossLimit) {
835 831
836 const G4double x = (fRange - length)/reduc << 832 G4double x = (fRange - length)/reduceFactor;
837 const G4double de = preStepKinEnergy - Sca << 833 //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl;
838 if(de > 0.0) { eloss = de; } << 834 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
>> 835
839 /* 836 /*
840 if(-1 < verboseLevel) 837 if(-1 < verboseLevel)
841 G4cout << " Long STEP: rPre(mm)=" << 838 G4cout << "Long STEP: rPre(mm)= "
842 << GetScaledRangeForScaledEnergy( 839 << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
843 << " x(mm)=" << x/mm << 840 << " rPost(mm)= " << x/mm
844 << " eloss(MeV)=" << eloss/MeV << 841 << " ePre(MeV)= " << preStepScaledEnergy/MeV
845 << " rFactor=" << reduceFactor << 842 << " eloss(MeV)= " << eloss/MeV << " eloss0(MeV)= "
846 << " massRatio=" << massRatio << 843 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
>> 844 << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
847 << G4endl; 845 << G4endl;
848 */ 846 */
849 } 847 }
850 848
851 /* 849 /*
852 if(-1 < verboseLevel ) { 850 if(-1 < verboseLevel ) {
853 G4cout << "Before fluct: eloss(MeV)= " << 851 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
854 << " e-eloss= " << preStepKinEnergy 852 << " e-eloss= " << preStepKinEnergy-eloss
855 << " step(mm)= " << length/mm << " 853 << " step(mm)= " << length/mm << " range(mm)= " << fRange/mm
856 << " fluct= " << lossFluctuationFla 854 << " fluct= " << lossFluctuationFlag << G4endl;
857 } 855 }
858 */ 856 */
859 857
860 const G4double cut = (*theCuts)[currentCoupl 858 const G4double cut = (*theCuts)[currentCoupleIndex];
861 G4double esec = 0.0; 859 G4double esec = 0.0;
862 860
863 // Corrections, which cannot be tabulated 861 // Corrections, which cannot be tabulated
864 if(isIon) { 862 if(isIon) {
865 currentModel->CorrectionsAlongStep(current 863 currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
866 length, 864 length, eloss);
867 eloss = std::max(eloss, 0.0); 865 eloss = std::max(eloss, 0.0);
868 } 866 }
869 867
870 // Sample fluctuations if not full energy lo 868 // Sample fluctuations if not full energy loss
871 if(eloss >= preStepKinEnergy) { 869 if(eloss >= preStepKinEnergy) {
872 eloss = preStepKinEnergy; 870 eloss = preStepKinEnergy;
873 871
874 } else if (lossFluctuationFlag) { 872 } else if (lossFluctuationFlag) {
875 const G4double tmax = currentModel->MaxSec 873 const G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
876 const G4double tcut = std::min(cut, tmax); 874 const G4double tcut = std::min(cut, tmax);
877 G4VEmFluctuationModel* fluc = currentModel 875 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
878 eloss = fluc->SampleFluctuations(currentCo 876 eloss = fluc->SampleFluctuations(currentCouple,dynParticle,
879 tcut, tma 877 tcut, tmax, length, eloss);
880 /* 878 /*
881 if(-1 < verboseLevel) 879 if(-1 < verboseLevel)
882 G4cout << "After fluct: eloss(MeV)= " << 880 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
883 << " fluc= " << (eloss-eloss0)/Me 881 << " fluc= " << (eloss-eloss0)/MeV
884 << " ChargeSqRatio= " << chargeSq 882 << " ChargeSqRatio= " << chargeSqRatio
885 << " massRatio= " << massRatio << 883 << " massRatio= " << massRatio << " tmax= " << tmax << G4endl;
886 */ 884 */
887 } 885 }
888 886
889 // deexcitation 887 // deexcitation
890 if (useDeexcitation) { 888 if (useDeexcitation) {
891 G4double esecfluo = preStepKinEnergy; 889 G4double esecfluo = preStepKinEnergy;
892 G4double de = esecfluo; 890 G4double de = esecfluo;
893 atomDeexcitation->AlongStepDeexcitation(sc 891 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
894 de 892 de, (G4int)currentCoupleIndex);
895 893
896 // sum of de-excitation energies 894 // sum of de-excitation energies
897 esecfluo -= de; 895 esecfluo -= de;
898 896
899 // subtracted from energy loss 897 // subtracted from energy loss
900 if(eloss >= esecfluo) { 898 if(eloss >= esecfluo) {
901 esec += esecfluo; 899 esec += esecfluo;
902 eloss -= esecfluo; 900 eloss -= esecfluo;
903 } else { 901 } else {
904 esec += esecfluo; 902 esec += esecfluo;
905 eloss = 0.0; 903 eloss = 0.0;
906 } 904 }
907 } 905 }
908 if(nullptr != subcutProducer && IsRegionForC 906 if(nullptr != subcutProducer && IsRegionForCubcutProcessor(track)) {
909 subcutProducer->SampleSecondaries(step, sc 907 subcutProducer->SampleSecondaries(step, scTracks, eloss, cut);
910 } 908 }
911 // secondaries from atomic de-excitation and 909 // secondaries from atomic de-excitation and subcut
912 if(!scTracks.empty()) { FillSecondariesAlong 910 if(!scTracks.empty()) { FillSecondariesAlongStep(weight); }
913 911
914 // Energy balance 912 // Energy balance
915 G4double finalT = preStepKinEnergy - eloss - 913 G4double finalT = preStepKinEnergy - eloss - esec;
916 if (finalT <= lowestKinEnergy) { 914 if (finalT <= lowestKinEnergy) {
917 eloss += finalT; 915 eloss += finalT;
918 finalT = 0.0; 916 finalT = 0.0;
919 } else if(isIon) { 917 } else if(isIon) {
920 fParticleChange.SetProposedCharge( 918 fParticleChange.SetProposedCharge(
921 currentModel->GetParticleCharge(track.Ge 919 currentModel->GetParticleCharge(track.GetParticleDefinition(),
922 currentM 920 currentMaterial,finalT));
923 } 921 }
924 eloss = std::max(eloss, 0.0); 922 eloss = std::max(eloss, 0.0);
925 923
926 fParticleChange.SetProposedKineticEnergy(fin 924 fParticleChange.SetProposedKineticEnergy(finalT);
927 fParticleChange.ProposeLocalEnergyDeposit(el 925 fParticleChange.ProposeLocalEnergyDeposit(eloss);
928 /* 926 /*
929 if(-1 < verboseLevel) { 927 if(-1 < verboseLevel) {
930 G4double del = finalT + eloss + esec - pre 928 G4double del = finalT + eloss + esec - preStepKinEnergy;
931 G4cout << "Final value eloss(MeV)= " << el 929 G4cout << "Final value eloss(MeV)= " << eloss/MeV
932 << " preStepKinEnergy= " << preStep 930 << " preStepKinEnergy= " << preStepKinEnergy
933 << " postStepKinEnergy= " << finalT 931 << " postStepKinEnergy= " << finalT
934 << " de(keV)= " << del/keV 932 << " de(keV)= " << del/keV
935 << " lossFlag= " << lossFluctuation 933 << " lossFlag= " << lossFluctuationFlag
936 << " status= " << track.GetTrackSt 934 << " status= " << track.GetTrackStatus()
937 << G4endl; 935 << G4endl;
938 } 936 }
939 */ 937 */
940 return &fParticleChange; 938 return &fParticleChange;
941 } 939 }
942 940
943 //....oooOO0OOooo........oooOO0OOooo........oo 941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
944 942
945 void G4VEnergyLossProcess::FillSecondariesAlon 943 void G4VEnergyLossProcess::FillSecondariesAlongStep(G4double wt)
946 { 944 {
947 const std::size_t n0 = scTracks.size(); 945 const std::size_t n0 = scTracks.size();
948 G4double weight = wt; 946 G4double weight = wt;
949 // weight may be changed by biasing manager 947 // weight may be changed by biasing manager
950 if(biasManager) { 948 if(biasManager) {
951 if(biasManager->SecondaryBiasingRegion((G4 949 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
952 weight *= 950 weight *=
953 biasManager->ApplySecondaryBiasing(scT 951 biasManager->ApplySecondaryBiasing(scTracks, (G4int)currentCoupleIndex);
954 } 952 }
955 } 953 }
956 954
957 // fill secondaries 955 // fill secondaries
958 const std::size_t n = scTracks.size(); 956 const std::size_t n = scTracks.size();
959 fParticleChange.SetNumberOfSecondaries((G4in 957 fParticleChange.SetNumberOfSecondaries((G4int)n);
960 958
961 for(std::size_t i=0; i<n; ++i) { 959 for(std::size_t i=0; i<n; ++i) {
962 G4Track* t = scTracks[i]; 960 G4Track* t = scTracks[i];
963 if(nullptr != t) { 961 if(nullptr != t) {
964 t->SetWeight(weight); 962 t->SetWeight(weight);
965 pParticleChange->AddSecondary(t); 963 pParticleChange->AddSecondary(t);
966 G4int pdg = t->GetDefinition()->GetPDGEn << 964 if(i >= n0) { t->SetCreatorModelID(biasID); }
967 if (i < n0) { <<
968 if (pdg == 22) { <<
969 t->SetCreatorModelID(gpixeID); <<
970 } else if (pdg == 11) { <<
971 t->SetCreatorModelID(epixeID); <<
972 } else { <<
973 t->SetCreatorModelID(biasID); <<
974 } <<
975 } else { <<
976 t->SetCreatorModelID(biasID); <<
977 } <<
978 } 965 }
979 } 966 }
980 scTracks.clear(); 967 scTracks.clear();
981 } 968 }
982 969
983 //....oooOO0OOooo........oooOO0OOooo........oo 970 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
984 971
985 G4VParticleChange* G4VEnergyLossProcess::PostS 972 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
986 973 const G4Step& step)
987 { 974 {
988 // clear number of interaction lengths in an 975 // clear number of interaction lengths in any case
989 theNumberOfInteractionLengthLeft = -1.0; 976 theNumberOfInteractionLengthLeft = -1.0;
990 mfpKinEnergy = DBL_MAX; 977 mfpKinEnergy = DBL_MAX;
991 978
992 fParticleChange.InitializeForPostStep(track) 979 fParticleChange.InitializeForPostStep(track);
993 const G4double finalT = track.GetKineticEner 980 const G4double finalT = track.GetKineticEnergy();
994 981
995 const G4double postStepScaledEnergy = finalT 982 const G4double postStepScaledEnergy = finalT*massRatio;
996 SelectModel(postStepScaledEnergy); 983 SelectModel(postStepScaledEnergy);
997 984
998 if(!currentModel->IsActive(postStepScaledEne 985 if(!currentModel->IsActive(postStepScaledEnergy)) {
999 return &fParticleChange; 986 return &fParticleChange;
1000 } 987 }
1001 /* 988 /*
1002 if(1 < verboseLevel) { 989 if(1 < verboseLevel) {
1003 G4cout<<GetProcessName()<<" PostStepDoIt: 990 G4cout<<GetProcessName()<<" PostStepDoIt: E(MeV)= "<< finalT/MeV<< G4endl;
1004 } 991 }
1005 */ 992 */
1006 // forced process - should happen only once 993 // forced process - should happen only once per track
1007 if(biasFlag) { 994 if(biasFlag) {
1008 if(biasManager->ForcedInteractionRegion(( 995 if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
1009 biasFlag = false; 996 biasFlag = false;
1010 } 997 }
1011 } 998 }
1012 const G4DynamicParticle* dp = track.GetDyna 999 const G4DynamicParticle* dp = track.GetDynamicParticle();
1013 1000
1014 // Integral approach 1001 // Integral approach
1015 if (fXSType != fEmNoIntegral) { 1002 if (fXSType != fEmNoIntegral) {
1016 const G4double logFinalT = dp->GetLogKine 1003 const G4double logFinalT = dp->GetLogKineticEnergy();
1017 G4double lx = GetLambdaForScaledEnergy(po 1004 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy,
1018 lo 1005 logFinalT + logMassRatio);
1019 lx = std::max(lx, 0.0); 1006 lx = std::max(lx, 0.0);
1020 1007
1021 // if both lg and lx are zero then no int 1008 // if both lg and lx are zero then no interaction
1022 if(preStepLambda*G4UniformRand() >= lx) { 1009 if(preStepLambda*G4UniformRand() >= lx) {
1023 return &fParticleChange; 1010 return &fParticleChange;
1024 } 1011 }
1025 } 1012 }
1026 1013
1027 // define new weight for primary and second 1014 // define new weight for primary and secondaries
1028 G4double weight = fParticleChange.GetParent 1015 G4double weight = fParticleChange.GetParentWeight();
1029 if(weightFlag) { 1016 if(weightFlag) {
1030 weight /= biasFactor; 1017 weight /= biasFactor;
1031 fParticleChange.ProposeWeight(weight); 1018 fParticleChange.ProposeWeight(weight);
1032 } 1019 }
1033 1020
1034 const G4double tcut = (*theCuts)[currentCou 1021 const G4double tcut = (*theCuts)[currentCoupleIndex];
1035 1022
1036 // sample secondaries 1023 // sample secondaries
1037 secParticles.clear(); 1024 secParticles.clear();
1038 currentModel->SampleSecondaries(&secParticl 1025 currentModel->SampleSecondaries(&secParticles, currentCouple, dp, tcut);
1039 1026
1040 const G4int num0 = (G4int)secParticles.size 1027 const G4int num0 = (G4int)secParticles.size();
1041 1028
1042 // bremsstrahlung splitting or Russian roul 1029 // bremsstrahlung splitting or Russian roulette
1043 if(biasManager) { 1030 if(biasManager) {
1044 if(biasManager->SecondaryBiasingRegion((G 1031 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
1045 G4double eloss = 0.0; 1032 G4double eloss = 0.0;
1046 weight *= biasManager->ApplySecondaryBi 1033 weight *= biasManager->ApplySecondaryBiasing(
1047 secPart 1034 secParticles,
1048 track, 1035 track, currentModel,
1049 &fParti 1036 &fParticleChange, eloss,
1050 (G4int) 1037 (G4int)currentCoupleIndex, tcut,
1051 step.Ge 1038 step.GetPostStepPoint()->GetSafety());
1052 if(eloss > 0.0) { 1039 if(eloss > 0.0) {
1053 eloss += fParticleChange.GetLocalEner 1040 eloss += fParticleChange.GetLocalEnergyDeposit();
1054 fParticleChange.ProposeLocalEnergyDep 1041 fParticleChange.ProposeLocalEnergyDeposit(eloss);
1055 } 1042 }
1056 } 1043 }
1057 } 1044 }
1058 1045
1059 // save secondaries 1046 // save secondaries
1060 const G4int num = (G4int)secParticles.size( 1047 const G4int num = (G4int)secParticles.size();
1061 if(num > 0) { 1048 if(num > 0) {
1062 1049
1063 fParticleChange.SetNumberOfSecondaries(nu 1050 fParticleChange.SetNumberOfSecondaries(num);
1064 G4double time = track.GetGlobalTime(); 1051 G4double time = track.GetGlobalTime();
1065 1052
1066 G4int n1(0), n2(0); 1053 G4int n1(0), n2(0);
1067 if(num0 > mainSecondaries) { 1054 if(num0 > mainSecondaries) {
1068 currentModel->FillNumberOfSecondaries(n 1055 currentModel->FillNumberOfSecondaries(n1, n2);
1069 } 1056 }
1070 1057
1071 for (G4int i=0; i<num; ++i) { 1058 for (G4int i=0; i<num; ++i) {
1072 if(nullptr != secParticles[i]) { 1059 if(nullptr != secParticles[i]) {
1073 G4Track* t = new G4Track(secParticles 1060 G4Track* t = new G4Track(secParticles[i], time, track.GetPosition());
1074 t->SetTouchableHandle(track.GetToucha 1061 t->SetTouchableHandle(track.GetTouchableHandle());
1075 if (biasManager) { 1062 if (biasManager) {
1076 t->SetWeight(weight * biasManager-> 1063 t->SetWeight(weight * biasManager->GetWeight(i));
1077 } else { 1064 } else {
1078 t->SetWeight(weight); 1065 t->SetWeight(weight);
1079 } 1066 }
1080 if(i < num0) { 1067 if(i < num0) {
1081 t->SetCreatorModelID(secID); 1068 t->SetCreatorModelID(secID);
1082 } else if(i < num0 + n1) { 1069 } else if(i < num0 + n1) {
1083 t->SetCreatorModelID(tripletID); 1070 t->SetCreatorModelID(tripletID);
1084 } else { 1071 } else {
1085 t->SetCreatorModelID(biasID); 1072 t->SetCreatorModelID(biasID);
1086 } 1073 }
1087 1074
1088 //G4cout << "Secondary(post step) has 1075 //G4cout << "Secondary(post step) has weight " << t->GetWeight()
1089 // << ", kenergy " << t->GetKin 1076 // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV"
1090 // << " time= " << time/ns << " 1077 // << " time= " << time/ns << " ns " << G4endl;
1091 pParticleChange->AddSecondary(t); 1078 pParticleChange->AddSecondary(t);
1092 } 1079 }
1093 } 1080 }
1094 } 1081 }
1095 1082
1096 if(0.0 == fParticleChange.GetProposedKineti 1083 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
1097 fAlive == fParticleChange.GetTrackStatus 1084 fAlive == fParticleChange.GetTrackStatus()) {
1098 if(particle->GetProcessManager()->GetAtRe 1085 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
1099 { fParticleChange.ProposeTrackStatus 1086 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
1100 else { fParticleChange.ProposeTrackStatus 1087 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
1101 } 1088 }
1102 1089
1103 /* 1090 /*
1104 if(-1 < verboseLevel) { 1091 if(-1 < verboseLevel) {
1105 G4cout << "::PostStepDoIt: Sample seconda 1092 G4cout << "::PostStepDoIt: Sample secondary; Efin= "
1106 << fParticleChange.GetProposedKineticEner 1093 << fParticleChange.GetProposedKineticEnergy()/MeV
1107 << " MeV; model= (" << currentMode 1094 << " MeV; model= (" << currentModel->LowEnergyLimit()
1108 << ", " << currentModel->HighEner 1095 << ", " << currentModel->HighEnergyLimit() << ")"
1109 << " preStepLambda= " << preStepL 1096 << " preStepLambda= " << preStepLambda
1110 << " dir= " << track.GetMomentumD 1097 << " dir= " << track.GetMomentumDirection()
1111 << " status= " << track.GetTrackS 1098 << " status= " << track.GetTrackStatus()
1112 << G4endl; 1099 << G4endl;
1113 } 1100 }
1114 */ 1101 */
1115 return &fParticleChange; 1102 return &fParticleChange;
1116 } 1103 }
1117 1104
1118 //....oooOO0OOooo........oooOO0OOooo........o 1105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1119 1106
1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl 1107 G4bool G4VEnergyLossProcess::StorePhysicsTable(
1121 const G4ParticleDefinition* part, cons 1108 const G4ParticleDefinition* part, const G4String& dir, G4bool ascii)
1122 { 1109 {
1123 if (!isMaster || nullptr != baseParticle || 1110 if (!isMaster || nullptr != baseParticle || part != particle ) return true;
1124 for(std::size_t i=0; i<7; ++i) { 1111 for(std::size_t i=0; i<7; ++i) {
1125 // ionisation table only for ionisation p << 1112 if(nullptr != theData->Table(i)) {
1126 if (nullptr == theData->Table(i) || (!isI << 1113 if(1 < verboseLevel) {
1127 continue; << 1114 G4cout << "G4VEnergyLossProcess::StorePhysicsTable i=" << i
1128 } << 1115 << " " << particle->GetParticleName()
1129 if (-1 < verboseLevel) { << 1116 << " " << GetProcessName()
1130 G4cout << "G4VEnergyLossProcess::StoreP << 1117 << " " << tnames[i] << " " << theData->Table(i) << G4endl;
1131 << " " << particle->GetParticleName() << 1118 }
1132 << " " << GetProcessName() << 1119 if(!G4EmTableUtil::StoreTable(this, part, theData->Table(i),
1133 << " " << tnames[i] << " " << theDat << 1120 dir, tnames[i], verboseLevel, ascii)) {
1134 } << 1121 return false;
1135 if (!G4EmTableUtil::StoreTable(this, part << 1122 }
1136 dir, tnames[i], verboseLevel, asci <<
1137 return false; <<
1138 } 1123 }
1139 } 1124 }
1140 return true; 1125 return true;
1141 } 1126 }
1142 1127
1143 //....oooOO0OOooo........oooOO0OOooo........o 1128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1144 1129
1145 G4bool 1130 G4bool
1146 G4VEnergyLossProcess::RetrievePhysicsTable(co 1131 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
1147 co 1132 const G4String& dir, G4bool ascii)
1148 { 1133 {
1149 if (!isMaster || nullptr != baseParticle || 1134 if (!isMaster || nullptr != baseParticle || part != particle ) return true;
1150 for(std::size_t i=0; i<7; ++i) { 1135 for(std::size_t i=0; i<7; ++i) {
1151 // ionisation table only for ionisation p <<
1152 if (!isIonisation && 1 == i) { continue; <<
1153 if(!G4EmTableUtil::RetrieveTable(this, pa 1136 if(!G4EmTableUtil::RetrieveTable(this, part, theData->Table(i), dir, tnames[i],
1154 verboseL 1137 verboseLevel, ascii, spline)) {
1155 return false; 1138 return false;
1156 } 1139 }
1157 } 1140 }
1158 return true; 1141 return true;
1159 } 1142 }
1160 1143
1161 //....oooOO0OOooo........oooOO0OOooo........o 1144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1162 1145
1163 G4double G4VEnergyLossProcess::GetDEDXDispers 1146 G4double G4VEnergyLossProcess::GetDEDXDispersion(
1164 const G4Mat 1147 const G4MaterialCutsCouple *couple,
1165 const G4Dyn 1148 const G4DynamicParticle* dp,
1166 G4dou 1149 G4double length)
1167 { 1150 {
1168 DefineMaterial(couple); 1151 DefineMaterial(couple);
1169 G4double ekin = dp->GetKineticEnergy(); 1152 G4double ekin = dp->GetKineticEnergy();
1170 SelectModel(ekin*massRatio); 1153 SelectModel(ekin*massRatio);
1171 G4double tmax = currentModel->MaxSecondaryK 1154 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
1172 G4double tcut = std::min(tmax,(*theCuts)[cu 1155 G4double tcut = std::min(tmax,(*theCuts)[currentCoupleIndex]);
1173 G4double d = 0.0; 1156 G4double d = 0.0;
1174 G4VEmFluctuationModel* fm = currentModel->G 1157 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
1175 if(nullptr != fm) { d = fm->Dispersion(curr 1158 if(nullptr != fm) { d = fm->Dispersion(currentMaterial,dp,tcut,tmax,length); }
1176 return d; 1159 return d;
1177 } 1160 }
1178 1161
1179 //....oooOO0OOooo........oooOO0OOooo........o 1162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1180 1163
1181 G4double 1164 G4double
1182 G4VEnergyLossProcess::CrossSectionPerVolume(G 1165 G4VEnergyLossProcess::CrossSectionPerVolume(G4double kineticEnergy,
1183 c 1166 const G4MaterialCutsCouple* couple,
1184 G 1167 G4double logKineticEnergy)
1185 { 1168 {
1186 // Cross section per volume is calculated 1169 // Cross section per volume is calculated
1187 DefineMaterial(couple); 1170 DefineMaterial(couple);
1188 G4double cross = 0.0; 1171 G4double cross = 0.0;
1189 if (nullptr != theLambdaTable) { 1172 if (nullptr != theLambdaTable) {
1190 cross = GetLambdaForScaledEnergy(kineticE 1173 cross = GetLambdaForScaledEnergy(kineticEnergy * massRatio,
1191 logKinet 1174 logKineticEnergy + logMassRatio);
1192 } else { 1175 } else {
1193 SelectModel(kineticEnergy*massRatio); 1176 SelectModel(kineticEnergy*massRatio);
1194 cross = (!baseMat) ? biasFactor : biasFac 1177 cross = (!baseMat) ? biasFactor : biasFactor*(*theDensityFactor)[currentCoupleIndex];
1195 cross *= (currentModel->CrossSectionPerVo 1178 cross *= (currentModel->CrossSectionPerVolume(currentMaterial, particle, kineticEnergy,
1196 1179 (*theCuts)[currentCoupleIndex]));
1197 } 1180 }
1198 return std::max(cross, 0.0); 1181 return std::max(cross, 0.0);
1199 } 1182 }
1200 1183
1201 //....oooOO0OOooo........oooOO0OOooo........o 1184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1202 1185
1203 G4double G4VEnergyLossProcess::MeanFreePath(c 1186 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
1204 { 1187 {
1205 DefineMaterial(track.GetMaterialCutsCouple( 1188 DefineMaterial(track.GetMaterialCutsCouple());
1206 const G4double kinEnergy = track.GetKine 1189 const G4double kinEnergy = track.GetKineticEnergy();
1207 const G4double logKinEnergy = track.GetDyna 1190 const G4double logKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
1208 const G4double cs = GetLambdaForScaledEnerg 1191 const G4double cs = GetLambdaForScaledEnergy(kinEnergy * massRatio,
1209 1192 logKinEnergy + logMassRatio);
1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1193 return (0.0 < cs) ? 1.0/cs : DBL_MAX;
1211 } 1194 }
1212 1195
1213 //....oooOO0OOooo........oooOO0OOooo........o 1196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1214 1197
1215 G4double G4VEnergyLossProcess::ContinuousStep 1198 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
1216 1199 G4double x, G4double y,
1217 1200 G4double& z)
1218 { 1201 {
1219 return AlongStepGetPhysicalInteractionLengt 1202 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &aGPILSelection);
1220 } 1203 }
1221 1204
1222 //....oooOO0OOooo........oooOO0OOooo........o 1205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1223 1206
1224 G4double G4VEnergyLossProcess::GetMeanFreePat 1207 G4double G4VEnergyLossProcess::GetMeanFreePath(
1225 const G4Track& t 1208 const G4Track& track,
1226 G4double, 1209 G4double,
1227 G4ForceCondition 1210 G4ForceCondition* condition)
1228 1211
1229 { 1212 {
1230 *condition = NotForced; 1213 *condition = NotForced;
1231 return MeanFreePath(track); 1214 return MeanFreePath(track);
1232 } 1215 }
1233 1216
1234 //....oooOO0OOooo........oooOO0OOooo........o 1217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1235 1218
1236 G4double G4VEnergyLossProcess::GetContinuousS 1219 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
1237 const G4Track&, 1220 const G4Track&,
1238 G4double, G4double, G4double& 1221 G4double, G4double, G4double&)
1239 { 1222 {
1240 return DBL_MAX; 1223 return DBL_MAX;
1241 } 1224 }
1242 1225
1243 //....oooOO0OOooo........oooOO0OOooo........o 1226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1244 1227
1245 G4PhysicsVector* 1228 G4PhysicsVector*
1246 G4VEnergyLossProcess::LambdaPhysicsVector(con 1229 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple,
1247 G4d 1230 G4double)
1248 { 1231 {
1249 DefineMaterial(couple); 1232 DefineMaterial(couple);
1250 G4PhysicsVector* v = (*theLambdaTable)[base 1233 G4PhysicsVector* v = (*theLambdaTable)[basedCoupleIndex];
1251 return new G4PhysicsVector(*v); 1234 return new G4PhysicsVector(*v);
1252 } 1235 }
1253 1236
1254 //....oooOO0OOooo........oooOO0OOooo........o 1237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1255 1238
1256 void 1239 void
1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT 1240 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
1258 { 1241 {
1259 if(1 < verboseLevel) { 1242 if(1 < verboseLevel) {
1260 G4cout << "### Set DEDX table " << p << " 1243 G4cout << "### Set DEDX table " << p << " " << theDEDXTable
1261 << " " << theDEDXunRestrictedTable << 1244 << " " << theDEDXunRestrictedTable << " " << theIonisationTable
1262 << " for " << particle->GetParticl 1245 << " for " << particle->GetParticleName()
1263 << " and process " << GetProcessNa 1246 << " and process " << GetProcessName()
1264 << " type=" << tType << " isIonisation:" 1247 << " type=" << tType << " isIonisation:" << isIonisation << G4endl;
1265 } 1248 }
1266 if(fTotal == tType) { 1249 if(fTotal == tType) {
1267 theDEDXunRestrictedTable = p; 1250 theDEDXunRestrictedTable = p;
1268 } else if(fRestricted == tType) { 1251 } else if(fRestricted == tType) {
1269 theDEDXTable = p; 1252 theDEDXTable = p;
1270 if(isMaster && nullptr == baseParticle) { 1253 if(isMaster && nullptr == baseParticle) {
1271 theData->UpdateTable(theDEDXTable, 0); 1254 theData->UpdateTable(theDEDXTable, 0);
1272 } 1255 }
1273 } else if(fIsIonisation == tType) { 1256 } else if(fIsIonisation == tType) {
1274 theIonisationTable = p; 1257 theIonisationTable = p;
1275 if(isMaster && nullptr == baseParticle) { 1258 if(isMaster && nullptr == baseParticle) {
1276 theData->UpdateTable(theIonisationTable 1259 theData->UpdateTable(theIonisationTable, 1);
1277 } 1260 }
1278 } 1261 }
1279 } 1262 }
1280 1263
1281 //....oooOO0OOooo........oooOO0OOooo........o 1264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1282 1265
1283 void G4VEnergyLossProcess::SetCSDARangeTable( 1266 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
1284 { 1267 {
1285 theCSDARangeTable = p; 1268 theCSDARangeTable = p;
1286 } 1269 }
1287 1270
1288 //....oooOO0OOooo........oooOO0OOooo........o 1271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1289 1272
1290 void G4VEnergyLossProcess::SetRangeTableForLo 1273 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
1291 { 1274 {
1292 theRangeTableForLoss = p; 1275 theRangeTableForLoss = p;
1293 } 1276 }
1294 1277
1295 //....oooOO0OOooo........oooOO0OOooo........o 1278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1296 1279
1297 void G4VEnergyLossProcess::SetInverseRangeTab 1280 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
1298 { 1281 {
1299 theInverseRangeTable = p; 1282 theInverseRangeTable = p;
1300 } 1283 }
1301 1284
1302 //....oooOO0OOooo........oooOO0OOooo........o 1285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1303 1286
1304 void G4VEnergyLossProcess::SetLambdaTable(G4P 1287 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
1305 { 1288 {
1306 if(1 < verboseLevel) { 1289 if(1 < verboseLevel) {
1307 G4cout << "### Set Lambda table " << p << 1290 G4cout << "### Set Lambda table " << p << " " << theLambdaTable
1308 << " for " << particle->GetParticl 1291 << " for " << particle->GetParticleName()
1309 << " and process " << GetProcessNa 1292 << " and process " << GetProcessName() << G4endl;
1310 } 1293 }
1311 theLambdaTable = p; 1294 theLambdaTable = p;
1312 tablesAreBuilt = true; 1295 tablesAreBuilt = true;
1313 1296
1314 if(isMaster && nullptr != p) { 1297 if(isMaster && nullptr != p) {
1315 delete theEnergyOfCrossSectionMax; 1298 delete theEnergyOfCrossSectionMax;
1316 theEnergyOfCrossSectionMax = nullptr; 1299 theEnergyOfCrossSectionMax = nullptr;
1317 if(fEmTwoPeaks == fXSType) { 1300 if(fEmTwoPeaks == fXSType) {
1318 if(nullptr != fXSpeaks) { 1301 if(nullptr != fXSpeaks) {
1319 for(auto & ptr : *fXSpeaks) { delete ptr; } 1302 for(auto & ptr : *fXSpeaks) { delete ptr; }
1320 delete fXSpeaks; 1303 delete fXSpeaks;
1321 } 1304 }
1322 G4LossTableBuilder* bld = lManager->Get 1305 G4LossTableBuilder* bld = lManager->GetTableBuilder();
1323 fXSpeaks = G4EmUtility::FillPeaksStruct 1306 fXSpeaks = G4EmUtility::FillPeaksStructure(p, bld);
1324 if(nullptr == fXSpeaks) { fXSType = fEm 1307 if(nullptr == fXSpeaks) { fXSType = fEmOnePeak; }
1325 } 1308 }
1326 if(fXSType == fEmOnePeak) { 1309 if(fXSType == fEmOnePeak) {
1327 theEnergyOfCrossSectionMax = G4EmUtilit 1310 theEnergyOfCrossSectionMax = G4EmUtility::FindCrossSectionMax(p);
1328 if(nullptr == theEnergyOfCrossSectionMa 1311 if(nullptr == theEnergyOfCrossSectionMax) { fXSType = fEmIncreasing; }
1329 } 1312 }
1330 } 1313 }
1331 } 1314 }
1332 1315
1333 //....oooOO0OOooo........oooOO0OOooo........o 1316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1334 1317
1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe 1318 void G4VEnergyLossProcess::SetEnergyOfCrossSectionMax(std::vector<G4double>* p)
1336 { 1319 {
1337 theEnergyOfCrossSectionMax = p; 1320 theEnergyOfCrossSectionMax = p;
1338 } 1321 }
1339 1322
1340 //....oooOO0OOooo........oooOO0OOooo........o 1323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1341 1324
1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: 1325 void G4VEnergyLossProcess::SetTwoPeaksXS(std::vector<G4TwoPeaksXS*>* ptr)
1343 { 1326 {
1344 fXSpeaks = ptr; 1327 fXSpeaks = ptr;
1345 } 1328 }
1346 1329
1347 //....oooOO0OOooo........oooOO0OOooo........o 1330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1348 1331
1349 const G4Element* G4VEnergyLossProcess::GetCur 1332 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
1350 { 1333 {
1351 return (nullptr != currentModel) 1334 return (nullptr != currentModel)
1352 ? currentModel->GetCurrentElement(current 1335 ? currentModel->GetCurrentElement(currentMaterial) : nullptr;
1353 } 1336 }
1354 1337
1355 //....oooOO0OOooo........oooOO0OOooo........o 1338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1356 1339
1357 void G4VEnergyLossProcess::SetCrossSectionBia 1340 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f,
1358 1341 G4bool flag)
1359 { 1342 {
1360 if(f > 0.0) { 1343 if(f > 0.0) {
1361 biasFactor = f; 1344 biasFactor = f;
1362 weightFlag = flag; 1345 weightFlag = flag;
1363 if(1 < verboseLevel) { 1346 if(1 < verboseLevel) {
1364 G4cout << "### SetCrossSectionBiasingFa 1347 G4cout << "### SetCrossSectionBiasingFactor: for "
1365 << " process " << GetProcessName 1348 << " process " << GetProcessName()
1366 << " biasFactor= " << f << " wei 1349 << " biasFactor= " << f << " weightFlag= " << flag
1367 << G4endl; 1350 << G4endl;
1368 } 1351 }
1369 } 1352 }
1370 } 1353 }
1371 1354
1372 //....oooOO0OOooo........oooOO0OOooo........o 1355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1373 1356
1374 void G4VEnergyLossProcess::ActivateForcedInte 1357 void G4VEnergyLossProcess::ActivateForcedInteraction(G4double length,
1375 1358 const G4String& region,
1376 1359 G4bool flag)
1377 { 1360 {
1378 if(nullptr == biasManager) { biasManager = 1361 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1379 if(1 < verboseLevel) { 1362 if(1 < verboseLevel) {
1380 G4cout << "### ActivateForcedInteraction: 1363 G4cout << "### ActivateForcedInteraction: for "
1381 << " process " << GetProcessName() 1364 << " process " << GetProcessName()
1382 << " length(mm)= " << length/mm 1365 << " length(mm)= " << length/mm
1383 << " in G4Region <" << region 1366 << " in G4Region <" << region
1384 << "> weightFlag= " << flag 1367 << "> weightFlag= " << flag
1385 << G4endl; 1368 << G4endl;
1386 } 1369 }
1387 weightFlag = flag; 1370 weightFlag = flag;
1388 biasManager->ActivateForcedInteraction(leng 1371 biasManager->ActivateForcedInteraction(length, region);
1389 } 1372 }
1390 1373
1391 //....oooOO0OOooo........oooOO0OOooo........o 1374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1392 1375
1393 void 1376 void
1394 G4VEnergyLossProcess::ActivateSecondaryBiasin 1377 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region,
1395 1378 G4double factor,
1396 1379 G4double energyLimit)
1397 { 1380 {
1398 if (0.0 <= factor) { 1381 if (0.0 <= factor) {
1399 // Range cut can be applied only for e- 1382 // Range cut can be applied only for e-
1400 if(0.0 == factor && secondaryParticle != 1383 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
1401 { return; } 1384 { return; }
1402 1385
1403 if(nullptr == biasManager) { biasManager 1386 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1404 biasManager->ActivateSecondaryBiasing(reg 1387 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
1405 if(1 < verboseLevel) { 1388 if(1 < verboseLevel) {
1406 G4cout << "### ActivateSecondaryBiasing 1389 G4cout << "### ActivateSecondaryBiasing: for "
1407 << " process " << GetProcessName 1390 << " process " << GetProcessName()
1408 << " factor= " << factor 1391 << " factor= " << factor
1409 << " in G4Region <" << region 1392 << " in G4Region <" << region
1410 << "> energyLimit(MeV)= " << ene 1393 << "> energyLimit(MeV)= " << energyLimit/MeV
1411 << G4endl; 1394 << G4endl;
1412 } 1395 }
1413 } 1396 }
1414 } 1397 }
1415 1398
1416 //....oooOO0OOooo........oooOO0OOooo........o 1399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1417 1400
1418 void G4VEnergyLossProcess::SetIonisation(G4bo 1401 void G4VEnergyLossProcess::SetIonisation(G4bool val)
1419 { 1402 {
1420 isIonisation = val; 1403 isIonisation = val;
1421 aGPILSelection = (val) ? CandidateForSelect 1404 aGPILSelection = (val) ? CandidateForSelection : NotCandidateForSelection;
1422 } 1405 }
1423 1406
1424 //....oooOO0OOooo........oooOO0OOooo........o 1407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1425 1408
1426 void G4VEnergyLossProcess::SetLinearLossLimi 1409 void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
1427 { 1410 {
1428 if(0.0 < val && val < 1.0) { 1411 if(0.0 < val && val < 1.0) {
1429 linLossLimit = val; 1412 linLossLimit = val;
1430 actLinLossLimit = true; 1413 actLinLossLimit = true;
1431 } else { PrintWarning("SetLinearLossLimit", 1414 } else { PrintWarning("SetLinearLossLimit", val); }
1432 } 1415 }
1433 1416
1434 //....oooOO0OOooo........oooOO0OOooo........o 1417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1435 1418
1436 void G4VEnergyLossProcess::SetStepFunction(G4 1419 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
1437 { 1420 {
1438 if(0.0 < v1 && 0.0 < v2) { 1421 if(0.0 < v1 && 0.0 < v2) {
1439 dRoverRange = std::min(1.0, v1); 1422 dRoverRange = std::min(1.0, v1);
1440 finalRange = std::min(v2, 1.e+50); 1423 finalRange = std::min(v2, 1.e+50);
1441 } else { 1424 } else {
1442 PrintWarning("SetStepFunctionV1", v1); 1425 PrintWarning("SetStepFunctionV1", v1);
1443 PrintWarning("SetStepFunctionV2", v2); 1426 PrintWarning("SetStepFunctionV2", v2);
1444 } 1427 }
1445 } 1428 }
1446 1429
1447 //....oooOO0OOooo........oooOO0OOooo........o 1430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1448 1431
1449 void G4VEnergyLossProcess::SetLowestEnergyLim 1432 void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val)
1450 { 1433 {
1451 if(1.e-18 < val && val < 1.e+50) { lowestKi 1434 if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; }
1452 else { PrintWarning("SetLowestEnergyLimit", 1435 else { PrintWarning("SetLowestEnergyLimit", val); }
1453 } 1436 }
1454 1437
1455 //....oooOO0OOooo........oooOO0OOooo........o 1438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1456 1439
1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i 1440 void G4VEnergyLossProcess::SetDEDXBinning(G4int n)
1458 { 1441 {
1459 if(2 < n && n < 1000000000) { 1442 if(2 < n && n < 1000000000) {
1460 nBins = n; 1443 nBins = n;
1461 actBinning = true; 1444 actBinning = true;
1462 } else { 1445 } else {
1463 G4double e = (G4double)n; 1446 G4double e = (G4double)n;
1464 PrintWarning("SetDEDXBinning", e); 1447 PrintWarning("SetDEDXBinning", e);
1465 } 1448 }
1466 } 1449 }
1467 1450
1468 //....oooOO0OOooo........oooOO0OOooo........o 1451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1469 1452
1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 1453 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
1471 { 1454 {
1472 if(1.e-18 < e && e < maxKinEnergy) { 1455 if(1.e-18 < e && e < maxKinEnergy) {
1473 minKinEnergy = e; 1456 minKinEnergy = e;
1474 actMinKinEnergy = true; 1457 actMinKinEnergy = true;
1475 } else { PrintWarning("SetMinKinEnergy", e) 1458 } else { PrintWarning("SetMinKinEnergy", e); }
1476 } 1459 }
1477 1460
1478 //....oooOO0OOooo........oooOO0OOooo........o 1461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1479 1462
1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 1463 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
1481 { 1464 {
1482 if(minKinEnergy < e && e < 1.e+50) { 1465 if(minKinEnergy < e && e < 1.e+50) {
1483 maxKinEnergy = e; 1466 maxKinEnergy = e;
1484 actMaxKinEnergy = true; 1467 actMaxKinEnergy = true;
1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS 1468 if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; }
1486 } else { PrintWarning("SetMaxKinEnergy", e) 1469 } else { PrintWarning("SetMaxKinEnergy", e); }
1487 } 1470 }
1488 1471
1489 //....oooOO0OOooo........oooOO0OOooo........o 1472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1490 1473
1491 void G4VEnergyLossProcess::PrintWarning(const 1474 void G4VEnergyLossProcess::PrintWarning(const G4String& tit, G4double val) const
1492 { 1475 {
1493 G4String ss = "G4VEnergyLossProcess::" + ti 1476 G4String ss = "G4VEnergyLossProcess::" + tit;
1494 G4ExceptionDescription ed; 1477 G4ExceptionDescription ed;
1495 ed << "Parameter is out of range: " << val 1478 ed << "Parameter is out of range: " << val
1496 << " it will have no effect!\n" << " Pr 1479 << " it will have no effect!\n" << " Process "
1497 << GetProcessName() << " nbins= " << nB 1480 << GetProcessName() << " nbins= " << nBins
1498 << " Emin(keV)= " << minKinEnergy/keV 1481 << " Emin(keV)= " << minKinEnergy/keV
1499 << " Emax(GeV)= " << maxKinEnergy/GeV; 1482 << " Emax(GeV)= " << maxKinEnergy/GeV;
1500 G4Exception(ss, "em0044", JustWarning, ed); 1483 G4Exception(ss, "em0044", JustWarning, ed);
1501 } 1484 }
1502 1485
1503 //....oooOO0OOooo........oooOO0OOooo........o 1486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1504 1487
1505 void G4VEnergyLossProcess::ProcessDescription 1488 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const
1506 { 1489 {
1507 if(nullptr != particle) { StreamInfo(out, * 1490 if(nullptr != particle) { StreamInfo(out, *particle, true); }
1508 } 1491 }
1509 1492
1510 //....oooOO0OOooo........oooOO0OOooo........o 1493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1511 1494