| 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" <<
63 #include "G4EmTableUtil.hh" <<
64 #include "G4VEmModel.hh" 62 #include "G4VEmModel.hh"
65 #include "G4VEmFluctuationModel.hh" 63 #include "G4VEmFluctuationModel.hh"
66 #include "G4DataVector.hh" 64 #include "G4DataVector.hh"
67 #include "G4PhysicsLogVector.hh" 65 #include "G4PhysicsLogVector.hh"
68 #include "G4VParticleChange.hh" 66 #include "G4VParticleChange.hh"
>> 67 #include "G4PhysicsModelCatalog.hh"
>> 68 #include "G4Gamma.hh"
69 #include "G4Electron.hh" 69 #include "G4Electron.hh"
>> 70 #include "G4Positron.hh"
70 #include "G4ProcessManager.hh" 71 #include "G4ProcessManager.hh"
71 #include "G4UnitsTable.hh" 72 #include "G4UnitsTable.hh"
>> 73 #include "G4ProductionCutsTable.hh"
72 #include "G4Region.hh" 74 #include "G4Region.hh"
73 #include "G4RegionStore.hh" 75 #include "G4RegionStore.hh"
74 #include "G4PhysicsTableHelper.hh" 76 #include "G4PhysicsTableHelper.hh"
75 #include "G4SafetyHelper.hh" 77 #include "G4SafetyHelper.hh"
76 #include "G4EmDataHandler.hh" <<
77 #include "G4TransportationManager.hh" 78 #include "G4TransportationManager.hh"
78 #include "G4VAtomDeexcitation.hh" 79 #include "G4VAtomDeexcitation.hh"
79 #include "G4VSubCutProducer.hh" 80 #include "G4VSubCutProducer.hh"
80 #include "G4EmBiasingManager.hh" 81 #include "G4EmBiasingManager.hh"
81 #include "G4Log.hh" 82 #include "G4Log.hh"
82 #include <iostream> 83 #include <iostream>
83 84
84 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 86
86 namespace <<
87 { <<
88 G4String tnames[7] = <<
89 {"DEDX","Ionisation","DEDXnr","CSDARange", <<
90 } <<
91 <<
92 <<
93 G4VEnergyLossProcess::G4VEnergyLossProcess(con 87 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name,
94 G4P 88 G4ProcessType type):
95 G4VContinuousDiscreteProcess(name, type) 89 G4VContinuousDiscreteProcess(name, type)
96 { 90 {
97 theParameters = G4EmParameters::Instance(); 91 theParameters = G4EmParameters::Instance();
98 SetVerboseLevel(1); 92 SetVerboseLevel(1);
99 93
100 // low energy limit 94 // low energy limit
101 lowestKinEnergy = theParameters->LowestElect 95 lowestKinEnergy = theParameters->LowestElectronEnergy();
102 96
103 // Size of tables 97 // Size of tables
104 minKinEnergy = 0.1*CLHEP::keV; 98 minKinEnergy = 0.1*CLHEP::keV;
105 maxKinEnergy = 100.0*CLHEP::TeV; 99 maxKinEnergy = 100.0*CLHEP::TeV;
106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; 100 maxKinEnergyCSDA = 1.0*CLHEP::GeV;
107 nBins = 84; 101 nBins = 84;
108 nBinsCSDA = 35; 102 nBinsCSDA = 35;
109 103
110 invLambdaFactor = 1.0/lambdaFactor; <<
111 <<
112 // default linear loss limit 104 // default linear loss limit
113 finalRange = 1.*CLHEP::mm; 105 finalRange = 1.*CLHEP::mm;
114 106
>> 107 // particle types
>> 108 theElectron = G4Electron::Electron();
>> 109 thePositron = G4Positron::Positron();
>> 110 theGamma = G4Gamma::Gamma();
>> 111
115 // run time objects 112 // run time objects
116 pParticleChange = &fParticleChange; 113 pParticleChange = &fParticleChange;
117 fParticleChange.SetSecondaryWeightByProcess( 114 fParticleChange.SetSecondaryWeightByProcess(true);
118 modelManager = new G4EmModelManager(); 115 modelManager = new G4EmModelManager();
119 safetyHelper = G4TransportationManager::GetT 116 safetyHelper = G4TransportationManager::GetTransportationManager()
120 ->GetSafetyHelper(); 117 ->GetSafetyHelper();
121 aGPILSelection = CandidateForSelection; 118 aGPILSelection = CandidateForSelection;
122 119
123 // initialise model 120 // initialise model
124 lManager = G4LossTableManager::Instance(); 121 lManager = G4LossTableManager::Instance();
125 lManager->Register(this); 122 lManager->Register(this);
126 isMaster = lManager->IsMaster(); <<
127 123
128 G4LossTableBuilder* bld = lManager->GetTable 124 G4LossTableBuilder* bld = lManager->GetTableBuilder();
129 theDensityFactor = bld->GetDensityFactors(); 125 theDensityFactor = bld->GetDensityFactors();
130 theDensityIdx = bld->GetCoupleIndexes(); 126 theDensityIdx = bld->GetCoupleIndexes();
131 127
132 scTracks.reserve(10); 128 scTracks.reserve(10);
133 secParticles.reserve(12); 129 secParticles.reserve(12);
134 emModels = new std::vector<G4VEmModel*>; <<
135 } 130 }
136 131
137 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 133
139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 134 G4VEnergyLossProcess::~G4VEnergyLossProcess()
140 { 135 {
141 if (isMaster) { << 136 /*
142 if(nullptr == baseParticle) { delete theDa << 137 G4cout << "** G4VEnergyLossProcess::~G4VEnergyLossProcess() for "
143 delete theEnergyOfCrossSectionMax; << 138 << GetProcessName() << " isMaster: " << isMaster
>> 139 << " basePart: " << baseParticle
>> 140 << G4endl;
>> 141 G4cout << " isIonisation " << isIonisation << " "
>> 142 << theDEDXTable << " " << theIonisationTable << G4endl;
>> 143 */
>> 144
>> 145 if (isMaster && nullptr == baseParticle) {
>> 146 if(nullptr != theDEDXTable) {
>> 147
>> 148 //G4cout << " theIonisationTable " << theIonisationTable << G4endl;
>> 149 if(theIonisationTable == theDEDXTable) { theIonisationTable = nullptr; }
>> 150 //G4cout << " delete theDEDXTable " << theDEDXTable << G4endl;
>> 151 theDEDXTable->clearAndDestroy();
>> 152 delete theDEDXTable;
>> 153 theDEDXTable = nullptr;
>> 154 }
>> 155 //G4cout << " theIonisationTable " << theIonisationTable << G4endl;
>> 156 if(nullptr != theIonisationTable) {
>> 157 //G4cout << " delete theIonisationTable " << theIonisationTable << G4endl;
>> 158 theIonisationTable->clearAndDestroy();
>> 159 delete theIonisationTable;
>> 160 theIonisationTable = nullptr;
>> 161 }
>> 162 if(nullptr != theDEDXunRestrictedTable && isIonisation) {
>> 163 theDEDXunRestrictedTable->clearAndDestroy();
>> 164 delete theDEDXunRestrictedTable;
>> 165 theDEDXunRestrictedTable = nullptr;
>> 166 }
>> 167 if(nullptr != theCSDARangeTable && isIonisation) {
>> 168 theCSDARangeTable->clearAndDestroy();
>> 169 delete theCSDARangeTable;
>> 170 theCSDARangeTable = nullptr;
>> 171 }
>> 172 //G4cout << "delete RangeTable: " << theRangeTableForLoss << G4endl;
>> 173 if(nullptr != theRangeTableForLoss && isIonisation) {
>> 174 theRangeTableForLoss->clearAndDestroy();
>> 175 delete theRangeTableForLoss;
>> 176 theRangeTableForLoss = nullptr;
>> 177 }
>> 178 //G4cout << "delete InvRangeTable: " << theInverseRangeTable << G4endl;
>> 179 if(nullptr != theInverseRangeTable && isIonisation /*&& !isIon*/) {
>> 180 theInverseRangeTable->clearAndDestroy();
>> 181 delete theInverseRangeTable;
>> 182 theInverseRangeTable = nullptr;
>> 183 }
>> 184 //G4cout << "delete LambdaTable: " << theLambdaTable << G4endl;
>> 185 if(nullptr != theLambdaTable) {
>> 186 theLambdaTable->clearAndDestroy();
>> 187 delete theLambdaTable;
>> 188 theLambdaTable = nullptr;
>> 189 }
144 if(nullptr != fXSpeaks) { 190 if(nullptr != fXSpeaks) {
145 for(auto const & v : *fXSpeaks) { delete 191 for(auto const & v : *fXSpeaks) { delete v; }
146 delete fXSpeaks; 192 delete fXSpeaks;
>> 193 fXSpeaks = nullptr;
147 } 194 }
148 } 195 }
>> 196 secParticles.clear();
149 delete modelManager; 197 delete modelManager;
150 delete biasManager; 198 delete biasManager;
151 delete scoffRegions; 199 delete scoffRegions;
152 delete emModels; 200 delete emModels;
153 lManager->DeRegister(this); 201 lManager->DeRegister(this);
>> 202 //G4cout << "** all removed" << G4endl;
154 } 203 }
155 204
156 //....oooOO0OOooo........oooOO0OOooo........oo 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 206
158 G4double G4VEnergyLossProcess::MinPrimaryEnerg 207 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
159 208 const G4Material*,
160 209 G4double cut)
161 { 210 {
162 return cut; 211 return cut;
163 } 212 }
164 213
165 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 215
167 void G4VEnergyLossProcess::AddEmModel(G4int or 216 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* ptr,
168 G4VEmFlu 217 G4VEmFluctuationModel* fluc,
169 const G4 218 const G4Region* region)
170 { 219 {
171 if(nullptr == ptr) { return; } 220 if(nullptr == ptr) { return; }
172 G4VEmFluctuationModel* afluc = (nullptr == f 221 G4VEmFluctuationModel* afluc = (nullptr == fluc) ? fluctModel : fluc;
173 modelManager->AddEmModel(order, ptr, afluc, 222 modelManager->AddEmModel(order, ptr, afluc, region);
174 ptr->SetParticleChange(pParticleChange, aflu 223 ptr->SetParticleChange(pParticleChange, afluc);
175 } 224 }
176 225
177 //....oooOO0OOooo........oooOO0OOooo........oo 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 227
179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod 228 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* ptr, G4int)
180 { 229 {
181 if(nullptr == ptr) { return; } 230 if(nullptr == ptr) { return; }
>> 231 if(nullptr == emModels) { emModels = new std::vector<G4VEmModel*>; }
182 if(!emModels->empty()) { 232 if(!emModels->empty()) {
183 for(auto & em : *emModels) { if(em == ptr) 233 for(auto & em : *emModels) { if(em == ptr) { return; } }
184 } 234 }
185 emModels->push_back(ptr); 235 emModels->push_back(ptr);
186 } 236 }
187 237
188 //....oooOO0OOooo........oooOO0OOooo........oo 238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 239
190 void G4VEnergyLossProcess::SetDynamicMassCharg 240 void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
191 241 G4double charge2ratio)
192 { 242 {
193 massRatio = massratio; 243 massRatio = massratio;
194 logMassRatio = G4Log(massRatio); 244 logMassRatio = G4Log(massRatio);
195 fFactor = charge2ratio*biasFactor; 245 fFactor = charge2ratio*biasFactor;
196 if(baseMat) { fFactor *= (*theDensityFactor) 246 if(baseMat) { fFactor *= (*theDensityFactor)[currentCoupleIndex]; }
197 chargeSqRatio = charge2ratio; 247 chargeSqRatio = charge2ratio;
198 reduceFactor = 1.0/(fFactor*massRatio); 248 reduceFactor = 1.0/(fFactor*massRatio);
199 } 249 }
200 250
201 //....oooOO0OOooo........oooOO0OOooo........oo 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 252
203 void 253 void
204 G4VEnergyLossProcess::PreparePhysicsTable(cons 254 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
205 { 255 {
206 particle = G4EmTableUtil::CheckIon(this, &pa << 256 if(1 < verboseLevel) {
207 verboseLe << 257 G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
>> 258 << GetProcessName() << " for " << part.GetParticleName()
>> 259 << " " << this << G4endl;
>> 260 }
>> 261 isMaster = lManager->IsMaster();
>> 262
>> 263 // Are particle defined?
>> 264 if(nullptr == particle) { particle = ∂ }
>> 265
>> 266 if(part.GetParticleType() == "nucleus") {
>> 267
>> 268 G4String pname = part.GetParticleName();
>> 269 if(pname != "deuteron" && pname != "triton" &&
>> 270 pname != "alpha+" && pname != "alpha") {
>> 271
>> 272 if(nullptr == theGenericIon) {
>> 273 theGenericIon =
>> 274 G4ParticleTable::GetParticleTable()->FindParticle("GenericIon");
>> 275 }
>> 276 isIon = true;
>> 277 if(particle != theGenericIon) {
>> 278 G4ProcessManager* pm = theGenericIon->GetProcessManager();
>> 279 G4ProcessVector* v = pm->GetAlongStepProcessVector();
>> 280 size_t n = v->size();
>> 281 for(size_t j=0; j<n; ++j) {
>> 282 if((*v)[j] == this) {
>> 283 particle = theGenericIon;
>> 284 break;
>> 285 }
>> 286 }
>> 287 }
>> 288 }
>> 289 }
208 290
209 if( particle != &part ) { 291 if( particle != &part ) {
210 if(!isIon) { lManager->RegisterExtraPartic << 292 if(!isIon) {
>> 293 lManager->RegisterExtraParticle(&part, this);
>> 294 }
211 if(1 < verboseLevel) { 295 if(1 < verboseLevel) {
212 G4cout << "### G4VEnergyLossProcess::Pre 296 G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()"
213 << " interrupted for " << GetProc << 297 << " interrupted for "
214 << part.GetParticleName() << " is << 298 << part.GetParticleName() << " isIon=" << isIon
215 << " spline=" << spline << G4endl << 299 << " baseMat=" << baseMat
>> 300 << " particle " << particle << " GenericIon " << theGenericIon
>> 301 << G4endl;
216 } 302 }
217 return; 303 return;
218 } 304 }
219 305
220 tablesAreBuilt = false; 306 tablesAreBuilt = false;
221 if (GetProcessSubType() == fIonisation) { Se <<
222 307
223 G4LossTableBuilder* bld = lManager->GetTable 308 G4LossTableBuilder* bld = lManager->GetTableBuilder();
224 lManager->PreparePhysicsTable(&part, this); << 309 lManager->PreparePhysicsTable(&part, this, isMaster);
225 310
226 // Base particle and set of models can be de 311 // Base particle and set of models can be defined here
227 InitialiseEnergyLossProcess(particle, basePa 312 InitialiseEnergyLossProcess(particle, baseParticle);
228 313
229 // parameters of the process 314 // parameters of the process
230 if(!actLossFluc) { lossFluctuationFlag = the 315 if(!actLossFluc) { lossFluctuationFlag = theParameters->LossFluctuation(); }
231 useCutAsFinalRange = theParameters->UseCutAs << 316 rndmStepFlag = theParameters->UseCutAsFinalRange();
232 if(!actMinKinEnergy) { minKinEnergy = thePar 317 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
233 if(!actMaxKinEnergy) { maxKinEnergy = thePar 318 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
234 if(!actBinning) { nBins = theParameters->Num 319 if(!actBinning) { nBins = theParameters->NumberOfBins(); }
235 maxKinEnergyCSDA = theParameters->MaxEnergyF 320 maxKinEnergyCSDA = theParameters->MaxEnergyForCSDARange();
236 nBinsCSDA = theParameters->NumberOfBinsPerDe 321 nBinsCSDA = theParameters->NumberOfBinsPerDecade()
237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi 322 *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy));
238 if(!actLinLossLimit) { linLossLimit = thePar 323 if(!actLinLossLimit) { linLossLimit = theParameters->LinearLossLimit(); }
239 lambdaFactor = theParameters->LambdaFactor() << 324 lambdaFactor = theParameters->LambdaFactor();
240 invLambdaFactor = 1.0/lambdaFactor; << 325 logLambdafactor = G4Log(lambdaFactor);
241 if(isMaster) { SetVerboseLevel(theParameters 326 if(isMaster) { SetVerboseLevel(theParameters->Verbose()); }
242 else { SetVerboseLevel(theParameters->Worker 327 else { SetVerboseLevel(theParameters->WorkerVerbose()); }
243 // integral option may be disabled <<
244 if(!theParameters->Integral()) { fXSType = f <<
245 328
246 theParameters->DefineRegParamForLoss(this); 329 theParameters->DefineRegParamForLoss(this);
247 330
248 fRangeEnergy = 0.0; << 331 fRangeEnergy = fLambdaEnergy = 0.0;
249 332
250 G4double initialCharge = particle->GetPDGCha 333 G4double initialCharge = particle->GetPDGCharge();
251 G4double initialMass = particle->GetPDGMas 334 G4double initialMass = particle->GetPDGMass();
252 335
253 theParameters->FillStepFunction(particle, th 336 theParameters->FillStepFunction(particle, this);
254 337
>> 338 // integral option may be disabled
>> 339 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; }
>> 340
255 // parameters for scaling from the base part 341 // parameters for scaling from the base particle
256 if (nullptr != baseParticle) { 342 if (nullptr != baseParticle) {
257 massRatio = (baseParticle->GetPDGMass() 343 massRatio = (baseParticle->GetPDGMass())/initialMass;
258 logMassRatio = G4Log(massRatio); 344 logMassRatio = G4Log(massRatio);
259 G4double q = initialCharge/baseParticle->G 345 G4double q = initialCharge/baseParticle->GetPDGCharge();
260 chargeSqRatio = q*q; 346 chargeSqRatio = q*q;
261 if(chargeSqRatio > 0.0) { reduceFactor = 1 347 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
262 } 348 }
263 lowestKinEnergy = (initialMass < CLHEP::MeV) 349 lowestKinEnergy = (initialMass < CLHEP::MeV)
264 ? theParameters->LowestElectronEnergy() 350 ? theParameters->LowestElectronEnergy()
265 : theParameters->LowestMuHadEnergy(); 351 : theParameters->LowestMuHadEnergy();
266 352
267 // Tables preparation 353 // Tables preparation
268 if (isMaster && nullptr == baseParticle) { 354 if (isMaster && nullptr == baseParticle) {
269 if(nullptr == theData) { theData = new G4E <<
270 355
271 if(nullptr != theDEDXTable && isIonisation 356 if(nullptr != theDEDXTable && isIonisation) {
272 if(nullptr != theIonisationTable && theD 357 if(nullptr != theIonisationTable && theDEDXTable != theIonisationTable) {
273 theData->CleanTable(0); << 358 theDEDXTable->clearAndDestroy();
274 theDEDXTable = theIonisationTable; << 359 delete theDEDXTable;
275 theIonisationTable = nullptr; << 360 theDEDXTable = theIonisationTable;
276 } << 361 }
277 } 362 }
278 363
279 theDEDXTable = theData->MakeTable(theDEDXT << 364 theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable);
280 bld->InitialiseBaseMaterials(theDEDXTable) 365 bld->InitialiseBaseMaterials(theDEDXTable);
281 theData->UpdateTable(theIonisationTable, 1 <<
282 366
283 if (theParameters->BuildCSDARange()) { 367 if (theParameters->BuildCSDARange()) {
284 theDEDXunRestrictedTable = theData->Make << 368 theDEDXunRestrictedTable =
285 if(isIonisation) { theCSDARangeTable = t << 369 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable);
>> 370 theCSDARangeTable =
>> 371 G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable);
286 } 372 }
287 373
288 theLambdaTable = theData->MakeTable(4); << 374 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
>> 375
289 if(isIonisation) { 376 if(isIonisation) {
290 theRangeTableForLoss = theData->MakeTabl << 377 theRangeTableForLoss =
291 theInverseRangeTable = theData->MakeTabl << 378 G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss);
>> 379 theInverseRangeTable =
>> 380 G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable);
292 } 381 }
293 } <<
294 382
>> 383 if(fXSType == fEmTwoPeaks) {
>> 384 const G4ProductionCutsTable* theCoupleTable=
>> 385 G4ProductionCutsTable::GetProductionCutsTable();
>> 386 size_t n = theCoupleTable->GetTableSize();
>> 387 if(nullptr == fXSpeaks) {
>> 388 fXSpeaks = new std::vector<G4TwoPeaksXS*>;
>> 389 }
>> 390 fXSpeaks->resize(n, nullptr);
>> 391 }
>> 392 }
>> 393 /*
>> 394 G4cout << "** G4VEnergyLossProcess::PreparePhysicsTable() for "
>> 395 << GetProcessName() << " and " << particle->GetParticleName()
>> 396 << " isMaster: " << isMaster << " isIonisation: "
>> 397 << isIonisation << G4endl;
>> 398 G4cout << " theDEDX: " << theDEDXTable
>> 399 << " theRange: " << theRangeTableForLoss
>> 400 << " theInverse: " << theInverseRangeTable
>> 401 << " theLambda: " << theLambdaTable << G4endl;
>> 402 */
295 // forced biasing 403 // forced biasing
296 if(nullptr != biasManager) { 404 if(nullptr != biasManager) {
297 biasManager->Initialise(part,GetProcessNam 405 biasManager->Initialise(part,GetProcessName(),verboseLevel);
298 biasFlag = false; 406 biasFlag = false;
299 } 407 }
>> 408
>> 409 // defined ID of secondary particles
>> 410 G4int stype = GetProcessSubType();
>> 411 if(stype == fBremsstrahlung) {
>> 412 secID = _Bremsstrahlung;
>> 413 biasID = _SplitBremsstrahlung;
>> 414 } else if(stype == fPairProdByCharged) {
>> 415 secID = _PairProduction;
>> 416 mainSecondaries = 2;
>> 417 }
300 baseMat = bld->GetBaseMaterialFlag(); 418 baseMat = bld->GetBaseMaterialFlag();
>> 419
>> 420 // initialisation of models
301 numberOfModels = modelManager->NumberOfModel 421 numberOfModels = modelManager->NumberOfModels();
302 currentModel = modelManager->GetModel(0); << 422 for(G4int i=0; i<numberOfModels; ++i) {
303 G4EmTableUtil::UpdateModels(this, modelManag << 423 G4VEmModel* mod = modelManager->GetModel(i);
304 numberOfModels, << 424 if(0 == i) { currentModel = mod; }
305 mainSecondaries, << 425 mod->SetMasterThread(isMaster);
306 theParameters->U << 426 mod->SetAngularGeneratorFlag(
307 theCuts = modelManager->Initialise(particle, << 427 theParameters->UseAngularGeneratorForIonisation());
308 verboseLe << 428 if(mod->HighEnergyLimit() > maxKinEnergy) {
>> 429 mod->SetHighEnergyLimit(maxKinEnergy);
>> 430 }
>> 431 mod->SetUseBaseMaterials(baseMat);
>> 432 SetEmModel(mod);
>> 433 }
>> 434 theCuts = modelManager->Initialise(particle, secondaryParticle,
>> 435 1.0, verboseLevel);
>> 436
309 // subcut processor 437 // subcut processor
310 if(isIonisation) { 438 if(isIonisation) {
311 subcutProducer = lManager->SubCutProducer( 439 subcutProducer = lManager->SubCutProducer();
312 } 440 }
313 if(1 == nSCoffRegions) { 441 if(1 == nSCoffRegions) {
314 if((*scoffRegions)[0]->GetName() == "Defau 442 if((*scoffRegions)[0]->GetName() == "DefaultRegionForTheWorld") {
315 delete scoffRegions; 443 delete scoffRegions;
316 scoffRegions = nullptr; 444 scoffRegions = nullptr;
317 nSCoffRegions = 0; 445 nSCoffRegions = 0;
318 } 446 }
319 } 447 }
320 448
321 if(1 < verboseLevel) { 449 if(1 < verboseLevel) {
322 G4cout << "G4VEnergyLossProcess::PrepearPh 450 G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done "
323 << " for " << GetProcessName() << " << 451 << " for local " << particle->GetParticleName()
324 << " isIon= " << isIon << " spline= << 452 << " isIon= " << isIon;
325 if(baseParticle) { 453 if(baseParticle) {
326 G4cout << "; base: " << baseParticle->Ge 454 G4cout << "; base: " << baseParticle->GetParticleName();
327 } 455 }
328 G4cout << G4endl; <<
329 G4cout << " chargeSqRatio= " << chargeSqRa 456 G4cout << " chargeSqRatio= " << chargeSqRatio
330 << " massRatio= " << massRatio 457 << " massRatio= " << massRatio
331 << " reduceFactor= " << reduceFacto 458 << " reduceFactor= " << reduceFactor << G4endl;
332 if (nSCoffRegions > 0) { 459 if (nSCoffRegions > 0) {
333 G4cout << " SubCut secondary production 460 G4cout << " SubCut secondary production is ON for regions: " << G4endl;
334 for (G4int i=0; i<nSCoffRegions; ++i) { 461 for (G4int i=0; i<nSCoffRegions; ++i) {
335 const G4Region* r = (*scoffRegions)[i] 462 const G4Region* r = (*scoffRegions)[i];
336 G4cout << " " << r->GetName( 463 G4cout << " " << r->GetName() << G4endl;
337 } 464 }
338 } else if(nullptr != subcutProducer) { 465 } else if(nullptr != subcutProducer) {
339 G4cout << " SubCut secondary production 466 G4cout << " SubCut secondary production is ON for all regions" << G4endl;
340 } 467 }
341 } 468 }
342 } 469 }
343 470
344 //....oooOO0OOooo........oooOO0OOooo........oo 471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 472
346 void G4VEnergyLossProcess::BuildPhysicsTable(c 473 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
347 { 474 {
348 if(1 < verboseLevel) { 475 if(1 < verboseLevel) {
349 G4cout << "### G4VEnergyLossProcess::Build 476 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
350 << GetProcessName() 477 << GetProcessName()
351 << " and particle " << part.GetPart 478 << " and particle " << part.GetParticleName()
352 << "; the first particle " << parti << 479 << "; local: " << particle->GetParticleName();
353 if(baseParticle) { 480 if(baseParticle) {
354 G4cout << "; base: " << baseParticle->Ge 481 G4cout << "; base: " << baseParticle->GetParticleName();
355 } 482 }
356 G4cout << G4endl; << 483 G4cout << " TablesAreBuilt= " << tablesAreBuilt
357 G4cout << " TablesAreBuilt= " << tables << 484 << " isIon= " << isIon << " " << this << G4endl;
358 << " spline=" << spline << " ptr: " <<
359 } 485 }
360 486
361 if(&part == particle) { 487 if(&part == particle) {
>> 488
362 if(isMaster) { 489 if(isMaster) {
363 lManager->BuildPhysicsTable(particle, th 490 lManager->BuildPhysicsTable(particle, this);
364 491
365 } else { 492 } else {
366 const auto masterProcess = << 493
>> 494 const G4VEnergyLossProcess* masterProcess =
367 static_cast<const G4VEnergyLossProcess 495 static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
368 496
369 numberOfModels = modelManager->NumberOfM << 497 // copy table pointers from master thread
370 G4EmTableUtil::BuildLocalElossProcess(th << 498 SetDEDXTable(masterProcess->DEDXTable(),fRestricted);
371 pa << 499 SetDEDXTable(masterProcess->DEDXunRestrictedTable(),fTotal);
372 tablesAreBuilt = true; << 500 SetDEDXTable(masterProcess->IonisationTable(),fIsIonisation);
>> 501 SetRangeTableForLoss(masterProcess->RangeTableForLoss());
>> 502 SetCSDARangeTable(masterProcess->CSDARangeTable());
>> 503 SetSecondaryRangeTable(masterProcess->SecondaryRangeTable());
>> 504 SetInverseRangeTable(masterProcess->InverseRangeTable());
>> 505 SetLambdaTable(masterProcess->LambdaTable());
>> 506 SetTwoPeaksXS(masterProcess->TwoPeaksXS());
>> 507 isIonisation = masterProcess->IsIonisationProcess();
373 baseMat = masterProcess->UseBaseMaterial 508 baseMat = masterProcess->UseBaseMaterial();
>> 509
>> 510 tablesAreBuilt = true;
>> 511 // local initialisation of models
>> 512 G4bool printing = true;
>> 513 for(G4int i=0; i<numberOfModels; ++i) {
>> 514 G4VEmModel* mod = GetModelByIndex(i, printing);
>> 515 G4VEmModel* mod0= masterProcess->GetModelByIndex(i, printing);
>> 516 mod->SetUseBaseMaterials(baseMat);
>> 517 mod->InitialiseLocal(particle, mod0);
>> 518 }
374 lManager->LocalPhysicsTables(particle, t 519 lManager->LocalPhysicsTables(particle, this);
375 } 520 }
376 521
377 // needs to be done only once 522 // needs to be done only once
378 safetyHelper->InitialiseHelper(); 523 safetyHelper->InitialiseHelper();
379 } 524 }
380 // Added tracking cut to avoid tracking arti 525 // Added tracking cut to avoid tracking artifacts
381 // and identified deexcitation flag 526 // and identified deexcitation flag
382 if(isIonisation) { 527 if(isIonisation) {
383 atomDeexcitation = lManager->AtomDeexcitat 528 atomDeexcitation = lManager->AtomDeexcitation();
384 if(nullptr != atomDeexcitation) { 529 if(nullptr != atomDeexcitation) {
385 if(atomDeexcitation->IsPIXEActive()) { u 530 if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; }
386 } 531 }
387 } 532 }
388 533
389 // protection against double printout 534 // protection against double printout
390 if(theParameters->IsPrintLocked()) { return; 535 if(theParameters->IsPrintLocked()) { return; }
391 536
392 // explicitly defined printout by particle n 537 // explicitly defined printout by particle name
393 G4String num = part.GetParticleName(); 538 G4String num = part.GetParticleName();
394 if(1 < verboseLevel || 539 if(1 < verboseLevel ||
395 (0 < verboseLevel && (num == "e-" || 540 (0 < verboseLevel && (num == "e-" ||
396 num == "e+" || n 541 num == "e+" || num == "mu+" ||
397 num == "mu-" || n 542 num == "mu-" || num == "proton"||
398 num == "pi+" || n 543 num == "pi+" || num == "pi-" ||
399 num == "kaon+" || n 544 num == "kaon+" || num == "kaon-" ||
400 num == "alpha" || n 545 num == "alpha" || num == "anti_proton" ||
401 num == "GenericIon" << 546 num == "GenericIon"|| num == "alpha+" )))
402 StreamInfo(G4cout, part); << 547 {
403 } << 548 StreamInfo(G4cout, part);
>> 549 }
>> 550
>> 551 /*
>> 552 G4cout << "** G4VEnergyLossProcess::BuildPhysicsTable() for "
>> 553 << GetProcessName() << " and " << particle->GetParticleName()
>> 554 << " isMaster: " << isMaster << " isIonisation: "
>> 555 << isIonisation << G4endl;
>> 556 G4cout << " theDEDX: " << theDEDXTable
>> 557 << " theRange: " << theRangeTableForLoss
>> 558 << " theInverse: " << theInverseRangeTable
>> 559 << " theLambda: " << theLambdaTable << G4endl;
>> 560 */
>> 561 //if(1 < verboseLevel || verb) {
404 if(1 < verboseLevel) { 562 if(1 < verboseLevel) {
405 G4cout << "### G4VEnergyLossProcess::Build 563 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
406 << GetProcessName() 564 << GetProcessName()
407 << " and particle " << part.GetPart 565 << " and particle " << part.GetParticleName();
408 if(isIonisation) { G4cout << " isIonisati 566 if(isIonisation) { G4cout << " isIonisation flag=1"; }
409 G4cout << " baseMat=" << baseMat << G4endl 567 G4cout << " baseMat=" << baseMat << G4endl;
410 } 568 }
411 } 569 }
412 570
413 //....oooOO0OOooo........oooOO0OOooo........oo 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414 572
415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED 573 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
416 { 574 {
>> 575 if(1 < verboseLevel ) {
>> 576 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
>> 577 << " for " << GetProcessName()
>> 578 << " and particle " << particle->GetParticleName()
>> 579 << G4endl;
>> 580 }
417 G4PhysicsTable* table = nullptr; 581 G4PhysicsTable* table = nullptr;
418 G4double emax = maxKinEnergy; 582 G4double emax = maxKinEnergy;
419 G4int bin = nBins; 583 G4int bin = nBins;
420 584
421 if(fTotal == tType) { 585 if(fTotal == tType) {
422 emax = maxKinEnergyCSDA; 586 emax = maxKinEnergyCSDA;
423 bin = nBinsCSDA; 587 bin = nBinsCSDA;
424 table = theDEDXunRestrictedTable; 588 table = theDEDXunRestrictedTable;
425 } else if(fRestricted == tType) { 589 } else if(fRestricted == tType) {
426 table = theDEDXTable; 590 table = theDEDXTable;
427 } else { 591 } else {
428 G4cout << "G4VEnergyLossProcess::BuildDEDX 592 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
429 << tType << G4endl; 593 << tType << G4endl;
430 } 594 }
>> 595
>> 596 // Access to materials
>> 597 const G4ProductionCutsTable* theCoupleTable=
>> 598 G4ProductionCutsTable::GetProductionCutsTable();
>> 599 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 600
431 if(1 < verboseLevel) { 601 if(1 < verboseLevel) {
432 G4cout << "G4VEnergyLossProcess::BuildDEDX << 602 G4cout << numOfCouples << " materials"
433 << " for " << GetProcessName() << 603 << " minKinEnergy= " << minKinEnergy
434 << " and " << particle->GetParticle << 604 << " maxKinEnergy= " << emax
435 << "spline=" << spline << G4endl; << 605 << " nbin= " << bin
>> 606 << " EmTableType= " << tType
>> 607 << " table= " << table << " " << this
>> 608 << G4endl;
436 } 609 }
437 if(nullptr == table) { return table; } 610 if(nullptr == table) { return table; }
438 611
439 G4LossTableBuilder* bld = lManager->GetTable 612 G4LossTableBuilder* bld = lManager->GetTableBuilder();
440 G4EmTableUtil::BuildDEDXTable(this, particle << 613 G4PhysicsLogVector* aVector = nullptr;
441 table, minKinE << 614 G4PhysicsLogVector* bVector = nullptr;
442 verboseLevel, << 615
>> 616 for(size_t i=0; i<numOfCouples; ++i) {
>> 617
>> 618 if(1 < verboseLevel) {
>> 619 G4cout << "G4VEnergyLossProcess::BuildDEDXVector Idx= " << i
>> 620 << " flagTable= " << table->GetFlag(i)
>> 621 << " Flag= " << bld->GetFlag(i) << G4endl;
>> 622 }
>> 623 if(bld->GetFlag(i)) {
>> 624
>> 625 // create physics vector and fill it
>> 626 const G4MaterialCutsCouple* couple =
>> 627 theCoupleTable->GetMaterialCutsCouple(i);
>> 628 if(nullptr != (*table)[i]) { delete (*table)[i]; }
>> 629 if(nullptr != bVector) {
>> 630 aVector = new G4PhysicsLogVector(*bVector);
>> 631 } else {
>> 632 bVector = new G4PhysicsLogVector(minKinEnergy, emax, bin, spline);
>> 633 aVector = bVector;
>> 634 }
>> 635
>> 636 modelManager->FillDEDXVector(aVector, couple, tType);
>> 637 if(spline) { aVector->FillSecondDerivatives(); }
>> 638
>> 639 // Insert vector for this material into the table
>> 640 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
>> 641 }
>> 642 }
>> 643
>> 644 if(1 < verboseLevel) {
>> 645 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
>> 646 << particle->GetParticleName()
>> 647 << " and process " << GetProcessName()
>> 648 << G4endl;
>> 649 if(2 < verboseLevel) G4cout << (*table) << G4endl;
>> 650 }
>> 651
443 return table; 652 return table;
444 } 653 }
445 654
446 //....oooOO0OOooo........oooOO0OOooo........oo 655 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 656
448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam << 657 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType)
449 { 658 {
450 if(nullptr == theLambdaTable) { return theLa << 659 G4PhysicsTable* table = nullptr;
451 660
452 G4double scale = theParameters->MaxKinEnergy << 661 if(fRestricted == tType) {
453 G4int nbin = << 662 table = theLambdaTable;
454 theParameters->NumberOfBinsPerDecade()*G4l << 663 } else {
455 scale = nbin/G4Log(scale); << 664 G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
456 << 665 << tType << G4endl;
>> 666 }
>> 667
>> 668 if(1 < verboseLevel) {
>> 669 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
>> 670 << tType << " for process "
>> 671 << GetProcessName() << " and particle "
>> 672 << particle->GetParticleName()
>> 673 << " EmTableType= " << tType
>> 674 << " table= " << table
>> 675 << G4endl;
>> 676 }
>> 677 if(nullptr == table) { return table; }
>> 678
>> 679 // Access to materials
>> 680 const G4ProductionCutsTable* theCoupleTable=
>> 681 G4ProductionCutsTable::GetProductionCutsTable();
>> 682 size_t numOfCouples = theCoupleTable->GetTableSize();
457 G4LossTableBuilder* bld = lManager->GetTable 683 G4LossTableBuilder* bld = lManager->GetTableBuilder();
458 G4EmTableUtil::BuildLambdaTable(this, partic << 684
459 bld, theLamb << 685 G4PhysicsLogVector* aVector = nullptr;
460 minKinEnergy << 686 G4double scale = G4Log(maxKinEnergy/minKinEnergy);
461 verboseLevel << 687
462 return theLambdaTable; << 688 for(size_t i=0; i<numOfCouples; ++i) {
>> 689
>> 690 if (bld->GetFlag(i)) {
>> 691
>> 692 // create physics vector and fill it
>> 693 const G4MaterialCutsCouple* couple =
>> 694 theCoupleTable->GetMaterialCutsCouple(i);
>> 695 delete (*table)[i];
>> 696
>> 697 G4bool startNull = true;
>> 698 G4double emin =
>> 699 MinPrimaryEnergy(particle,couple->GetMaterial(),(*theCuts)[i]);
>> 700 if(minKinEnergy > emin) {
>> 701 emin = minKinEnergy;
>> 702 startNull = false;
>> 703 }
>> 704
>> 705 G4double emax = maxKinEnergy;
>> 706 if(emax <= emin) { emax = 2*emin; }
>> 707 G4int bin = G4lrint(nBins*G4Log(emax/emin)/scale);
>> 708 bin = std::max(bin, 3);
>> 709 aVector = new G4PhysicsLogVector(emin, emax, bin, spline);
>> 710
>> 711 modelManager->FillLambdaVector(aVector, couple, startNull, tType);
>> 712 if(spline) { aVector->FillSecondDerivatives(); }
>> 713
>> 714 // Insert vector for this material into the table
>> 715 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
>> 716 }
>> 717 }
>> 718
>> 719 if(1 < verboseLevel) {
>> 720 G4cout << "Lambda table is built for "
>> 721 << particle->GetParticleName()
>> 722 << G4endl;
>> 723 }
>> 724
>> 725 return table;
463 } 726 }
464 727
465 //....oooOO0OOooo........oooOO0OOooo........oo 728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 729
467 void G4VEnergyLossProcess::StreamInfo(std::ost 730 void G4VEnergyLossProcess::StreamInfo(std::ostream& out,
468 const G4ParticleDefinition& pa 731 const G4ParticleDefinition& part, G4bool rst) const
469 { 732 {
470 G4String indent = (rst ? " " : ""); 733 G4String indent = (rst ? " " : "");
471 out << std::setprecision(6); 734 out << std::setprecision(6);
472 out << G4endl << indent << GetProcessName() 735 out << G4endl << indent << GetProcessName() << ": ";
473 if (!rst) out << " for " << part.GetParticle 736 if (!rst) out << " for " << part.GetParticleName();
474 out << " XStype:" << fXSType 737 out << " XStype:" << fXSType
475 << " SubType=" << GetProcessSubType() < 738 << " SubType=" << GetProcessSubType() << G4endl
476 << " dE/dx and range tables from " 739 << " dE/dx and range tables from "
477 << G4BestUnit(minKinEnergy,"Energy") 740 << G4BestUnit(minKinEnergy,"Energy")
478 << " to " << G4BestUnit(maxKinEnergy,"En 741 << " to " << G4BestUnit(maxKinEnergy,"Energy")
479 << " in " << nBins << " bins" << G4endl 742 << " in " << nBins << " bins" << G4endl
480 << " Lambda tables from threshold t 743 << " Lambda tables from threshold to "
481 << G4BestUnit(maxKinEnergy,"Energy") 744 << G4BestUnit(maxKinEnergy,"Energy")
482 << ", " << theParameters->NumberOfBinsPe 745 << ", " << theParameters->NumberOfBinsPerDecade()
483 << " bins/decade, spline: " << spline 746 << " bins/decade, spline: " << spline
484 << G4endl; 747 << G4endl;
485 if(nullptr != theRangeTableForLoss && isIoni 748 if(nullptr != theRangeTableForLoss && isIonisation) {
486 out << " StepFunction=(" << dRoverRan 749 out << " StepFunction=(" << dRoverRange << ", "
487 << finalRange/mm << " mm)" 750 << finalRange/mm << " mm)"
488 << ", integ: " << fXSType 751 << ", integ: " << fXSType
489 << ", fluct: " << lossFluctuationFlag 752 << ", fluct: " << lossFluctuationFlag
490 << ", linLossLim= " << linLossLimit 753 << ", linLossLim= " << linLossLimit
491 << G4endl; 754 << G4endl;
492 } 755 }
493 StreamProcessInfo(out); 756 StreamProcessInfo(out);
494 modelManager->DumpModelList(out, verboseLeve 757 modelManager->DumpModelList(out, verboseLevel);
495 if(nullptr != theCSDARangeTable && isIonisat 758 if(nullptr != theCSDARangeTable && isIonisation) {
496 out << " CSDA range table up" 759 out << " CSDA range table up"
497 << " to " << G4BestUnit(maxKinEnergyCS 760 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
498 << " in " << nBinsCSDA << " bins" << G 761 << " in " << nBinsCSDA << " bins" << G4endl;
499 } 762 }
500 if(nSCoffRegions>0 && isIonisation) { 763 if(nSCoffRegions>0 && isIonisation) {
501 out << " Subcutoff sampling in " << n 764 out << " Subcutoff sampling in " << nSCoffRegions
502 << " regions" << G4endl; 765 << " regions" << G4endl;
503 } 766 }
504 if(2 < verboseLevel) { 767 if(2 < verboseLevel) {
505 for(std::size_t i=0; i<7; ++i) { << 768 out << " DEDXTable address= " << theDEDXTable << G4endl;
506 auto ta = theData->Table(i); << 769 if(nullptr != theDEDXTable && isIonisation)
507 out << " " << tnames[i] << " addres << 770 out << (*theDEDXTable) << G4endl;
508 if(nullptr != ta) { out << *ta << G4endl << 771 out << "non restricted DEDXTable address= "
>> 772 << theDEDXunRestrictedTable << G4endl;
>> 773 if(nullptr != theDEDXunRestrictedTable && isIonisation) {
>> 774 out << (*theDEDXunRestrictedTable) << G4endl;
>> 775 }
>> 776 out << " CSDARangeTable address= " << theCSDARangeTable << G4endl;
>> 777 if(nullptr != theCSDARangeTable && isIonisation) {
>> 778 out << (*theCSDARangeTable) << G4endl;
>> 779 }
>> 780 out << " RangeTableForLoss address= " << theRangeTableForLoss
>> 781 << G4endl;
>> 782 if(nullptr != theRangeTableForLoss && isIonisation) {
>> 783 out << (*theRangeTableForLoss) << G4endl;
>> 784 }
>> 785 out << " InverseRangeTable address= " << theInverseRangeTable
>> 786 << G4endl;
>> 787 if(nullptr != theInverseRangeTable && isIonisation) {
>> 788 out << (*theInverseRangeTable) << G4endl;
>> 789 }
>> 790 out << " LambdaTable address= " << theLambdaTable << G4endl;
>> 791 if(nullptr != theLambdaTable) {
>> 792 out << (*theLambdaTable) << G4endl;
509 } 793 }
510 } 794 }
511 } 795 }
512 796
513 //....oooOO0OOooo........oooOO0OOooo........oo 797 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
514 798
515 void G4VEnergyLossProcess::ActivateSubCutoff(c 799 void G4VEnergyLossProcess::ActivateSubCutoff(const G4Region* r)
516 { 800 {
517 if(nullptr == scoffRegions) { 801 if(nullptr == scoffRegions) {
518 scoffRegions = new std::vector<const G4Reg 802 scoffRegions = new std::vector<const G4Region*>;
519 } 803 }
520 // the region is in the list 804 // the region is in the list
521 if(!scoffRegions->empty()) { 805 if(!scoffRegions->empty()) {
522 for (auto & reg : *scoffRegions) { 806 for (auto & reg : *scoffRegions) {
523 if (reg == r) { return; } 807 if (reg == r) { return; }
524 } 808 }
525 } 809 }
526 // new region 810 // new region
527 scoffRegions->push_back(r); 811 scoffRegions->push_back(r);
528 ++nSCoffRegions; 812 ++nSCoffRegions;
529 } 813 }
530 814
531 //....oooOO0OOooo........oooOO0OOooo........oo 815 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
532 816
533 G4bool G4VEnergyLossProcess::IsRegionForCubcut 817 G4bool G4VEnergyLossProcess::IsRegionForCubcutProcessor(const G4Track& aTrack)
534 { 818 {
535 if(0 == nSCoffRegions) { return true; } 819 if(0 == nSCoffRegions) { return true; }
536 const G4Region* r = aTrack.GetVolume()->GetL 820 const G4Region* r = aTrack.GetVolume()->GetLogicalVolume()->GetRegion();
537 for(auto & reg : *scoffRegions) { 821 for(auto & reg : *scoffRegions) {
538 if(r == reg) { return true; } 822 if(r == reg) { return true; }
539 } 823 }
540 return false; 824 return false;
541 } 825 }
542 826
543 //....oooOO0OOooo........oooOO0OOooo........oo 827 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544 828
545 void G4VEnergyLossProcess::StartTracking(G4Tra 829 void G4VEnergyLossProcess::StartTracking(G4Track* track)
546 { 830 {
>> 831 /*
>> 832 G4cout << "G4VEnergyLossProcess::StartTracking: "
>> 833 << track->GetDefinition()->GetParticleName()
>> 834 << " e(MeV)= " << track->GetKineticEnergy();
>> 835 if(particle) G4cout << " " << particle->GetParticleName();
>> 836 if(baseParticle) G4cout << " basePart: " << baseParticle->GetParticleName();
>> 837 G4cout << " " << GetProcessName();
>> 838 if(isIon) G4cout << " isIon: Q=" << track->GetDefinition()->GetPDGCharge()
>> 839 << " Qdyn=" << track->GetDynamicParticle()->GetCharge();
>> 840 G4cout << G4endl;
>> 841 */
547 // reset parameters for the new track 842 // reset parameters for the new track
548 theNumberOfInteractionLengthLeft = -1.0; 843 theNumberOfInteractionLengthLeft = -1.0;
549 mfpKinEnergy = DBL_MAX; 844 mfpKinEnergy = DBL_MAX;
550 preStepLambda = 0.0; <<
551 currentCouple = nullptr; 845 currentCouple = nullptr;
552 846
553 // reset ion 847 // reset ion
554 if(isIon) { 848 if(isIon) {
555 const G4double newmass = track->GetDefinit 849 const G4double newmass = track->GetDefinition()->GetPDGMass();
556 massRatio = (nullptr == baseParticle) ? CL << 850 if(nullptr != baseParticle) {
557 : baseParticle->GetPDGMass()/newmass; << 851 massRatio = baseParticle->GetPDGMass()/newmass;
558 logMassRatio = G4Log(massRatio); << 852 logMassRatio = G4Log(massRatio);
>> 853 } else if(nullptr != theGenericIon) {
>> 854 massRatio = CLHEP::proton_mass_c2/newmass;
>> 855 logMassRatio = G4Log(massRatio);
>> 856 } else {
>> 857 massRatio = 1.0;
>> 858 logMassRatio = 0.0;
>> 859 }
559 } 860 }
560 // forced biasing only for primary particles 861 // forced biasing only for primary particles
561 if(nullptr != biasManager) { << 862 if(biasManager) {
562 if(0 == track->GetParentID()) { 863 if(0 == track->GetParentID()) {
563 biasFlag = true; 864 biasFlag = true;
564 biasManager->ResetForcedInteraction(); 865 biasManager->ResetForcedInteraction();
565 } 866 }
566 } 867 }
567 } 868 }
568 869
569 //....oooOO0OOooo........oooOO0OOooo........oo 870 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
570 871
571 G4double G4VEnergyLossProcess::AlongStepGetPhy 872 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
572 const G4Track& tr << 873 const G4Track&,G4double,G4double,G4double&,
573 G4GPILSelection* 874 G4GPILSelection* selection)
574 { 875 {
575 G4double x = DBL_MAX; 876 G4double x = DBL_MAX;
576 *selection = aGPILSelection; 877 *selection = aGPILSelection;
577 if(isIonisation && currentModel->IsActive(pr 878 if(isIonisation && currentModel->IsActive(preStepScaledEnergy)) {
578 GetScaledRangeForScaledEnergy(preStepScale << 879 GetScaledRangeForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy);
579 x = (useCutAsFinalRange) ? std::min(finalR << 880 const G4double finR = (rndmStepFlag) ? std::min(finalRange,
580 currentCouple->GetProductionCuts()->GetP 881 currentCouple->GetProductionCuts()->GetProductionCut(1)) : finalRange;
581 x = (fRange > x) ? fRange*dRoverRange + x* << 882 x = (fRange > finR) ?
582 : fRange; << 883 fRange*dRoverRange + finR*(1.0-dRoverRange)*(2.0-finR/fRange) : fRange;
583 /* << 884 // if(particle->GetPDGMass() > 0.9*GeV)
584 G4cout<<"AlongStepGPIL: " << GetProcessN << 885 /*
585 << " fRange=" << fRange << " finR=" << finR << 886 G4cout << GetProcessName() << ": e= " << preStepKinEnergy
>> 887 <<" range= "<<fRange << " idx= " << basedCoupleIndex
>> 888 << " finR= " << finR << " limit= " << x <<
>> 889 << "\n" << "massRatio= " << massRatio << " Q^2= " << chargeSqRatio
>> 890 << " dRoverRange= " << dRoverRange
>> 891 << " finalRange= " << finalRange << G4endl;
586 */ 892 */
587 } 893 }
>> 894 //G4cout<<"AlongStepGPIL: " << GetProcessName()<<": e= "<<preStepKinEnergy
>> 895 //<<" stepLimit= "<<x<<G4endl;
588 return x; 896 return x;
589 } 897 }
590 898
591 //....oooOO0OOooo........oooOO0OOooo........oo 899 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
592 900
593 G4double G4VEnergyLossProcess::PostStepGetPhys 901 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
594 const G4Track& tr 902 const G4Track& track,
595 G4double previo 903 G4double previousStepSize,
596 G4ForceCondition* 904 G4ForceCondition* condition)
597 { 905 {
598 // condition is set to "Not Forced" 906 // condition is set to "Not Forced"
599 *condition = NotForced; 907 *condition = NotForced;
600 G4double x = DBL_MAX; 908 G4double x = DBL_MAX;
601 909
602 // initialisation of material, mass, charge, 910 // initialisation of material, mass, charge, model
603 // at the beginning of the step 911 // at the beginning of the step
604 DefineMaterial(track.GetMaterialCutsCouple() 912 DefineMaterial(track.GetMaterialCutsCouple());
605 preStepKinEnergy = track.GetKineticEne 913 preStepKinEnergy = track.GetKineticEnergy();
>> 914 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
606 preStepScaledEnergy = preStepKinEnergy*ma 915 preStepScaledEnergy = preStepKinEnergy*massRatio;
>> 916 preStepLogScaledEnergy = preStepLogKinEnergy + logMassRatio;
607 SelectModel(preStepScaledEnergy); 917 SelectModel(preStepScaledEnergy);
608 918
609 if(!currentModel->IsActive(preStepScaledEner 919 if(!currentModel->IsActive(preStepScaledEnergy)) {
610 theNumberOfInteractionLengthLeft = -1.0; 920 theNumberOfInteractionLengthLeft = -1.0;
611 mfpKinEnergy = DBL_MAX; <<
612 preStepLambda = 0.0; <<
613 currentInteractionLength = DBL_MAX; 921 currentInteractionLength = DBL_MAX;
614 return x; << 922 return x;
615 } 923 }
616 924
617 // change effective charge of a charged part 925 // change effective charge of a charged particle on fly
618 if(isIon) { 926 if(isIon) {
619 const G4double q2 = currentModel->ChargeSq 927 const G4double q2 = currentModel->ChargeSquareRatio(track);
620 fFactor = q2*biasFactor; << 928 if(q2 != chargeSqRatio) {
621 if(baseMat) { fFactor *= (*theDensityFacto << 929 fFactor *= q2/chargeSqRatio;
622 reduceFactor = 1.0/(fFactor*massRatio); << 930 reduceFactor = 1.0/(fFactor*massRatio);
623 if (lossFluctuationFlag) { << 931 chargeSqRatio = q2;
624 auto fluc = currentModel->GetModelOfFluc << 932 // G4cout << "PostStepGPIL: Q^2=" << chargeSqRatio << " reducedFactor=" << reduceFactor << G4endl;
625 fluc->SetParticleAndCharge(track.GetDefi <<
626 } 933 }
627 } 934 }
628 935
629 // forced biasing only for primary particles 936 // forced biasing only for primary particles
630 if(biasManager) { 937 if(biasManager) {
631 if(0 == track.GetParentID() && biasFlag && 938 if(0 == track.GetParentID() && biasFlag &&
632 biasManager->ForcedInteractionRegion((G << 939 biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
633 return biasManager->GetStepLimit((G4int) << 940 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
634 } 941 }
635 } 942 }
636 943
637 ComputeLambdaForScaledEnergy(preStepScaledEn << 944 // compute mean free path
638 << 945 ComputeLambdaForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy);
>> 946
639 // zero cross section 947 // zero cross section
640 if(preStepLambda <= 0.0) { 948 if(preStepLambda <= 0.0) {
641 theNumberOfInteractionLengthLeft = -1.0; 949 theNumberOfInteractionLengthLeft = -1.0;
642 currentInteractionLength = DBL_MAX; 950 currentInteractionLength = DBL_MAX;
643 } else { 951 } else {
644 952
645 // non-zero cross section 953 // non-zero cross section
646 if (theNumberOfInteractionLengthLeft < 0.0 954 if (theNumberOfInteractionLengthLeft < 0.0) {
647 955
648 // beggining of tracking (or just after 956 // beggining of tracking (or just after DoIt of this process)
649 theNumberOfInteractionLengthLeft = -G4Lo 957 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
650 theInitialNumberOfInteractionLength = th 958 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
651 959
652 } else if(currentInteractionLength < DBL_M 960 } else if(currentInteractionLength < DBL_MAX) {
653 961
654 // subtract NumberOfInteractionLengthLef 962 // subtract NumberOfInteractionLengthLeft using previous step
655 theNumberOfInteractionLengthLeft -= 963 theNumberOfInteractionLengthLeft -=
656 previousStepSize/currentInteractionLen 964 previousStepSize/currentInteractionLength;
657 965
658 theNumberOfInteractionLengthLeft = 966 theNumberOfInteractionLengthLeft =
659 std::max(theNumberOfInteractionLengthL 967 std::max(theNumberOfInteractionLengthLeft, 0.0);
660 } 968 }
661 969
662 // new mean free path and step limit 970 // new mean free path and step limit
663 currentInteractionLength = 1.0/preStepLamb 971 currentInteractionLength = 1.0/preStepLambda;
664 x = theNumberOfInteractionLengthLeft * cur 972 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
665 } 973 }
666 #ifdef G4VERBOSE 974 #ifdef G4VERBOSE
667 if (verboseLevel>2) { << 975 if (verboseLevel>2){
>> 976 // if(particle->GetPDGMass() > 0.9*GeV){
668 G4cout << "G4VEnergyLossProcess::PostStepG 977 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
669 G4cout << "[ " << GetProcessName() << "]" 978 G4cout << "[ " << GetProcessName() << "]" << G4endl;
670 G4cout << " for " << track.GetDefinition() 979 G4cout << " for " << track.GetDefinition()->GetParticleName()
671 << " in Material " << currentMate 980 << " in Material " << currentMaterial->GetName()
672 << " Ekin(MeV)= " << preStepKinEner 981 << " Ekin(MeV)= " << preStepKinEnergy/MeV
673 << " track material: " << track.Get << 982 << " " << track.GetMaterial()->GetName()
674 <<G4endl; 983 <<G4endl;
675 G4cout << "MeanFreePath = " << currentInte 984 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
676 << "InteractionLength= " << x/cm << 985 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
677 } 986 }
678 #endif 987 #endif
679 return x; 988 return x;
680 } 989 }
681 990
682 //....oooOO0OOooo........oooOO0OOooo........oo 991 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683 992
684 void 993 void
685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 994 G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e, G4double loge)
686 { 995 {
687 // cross section increased with energy 996 // cross section increased with energy
688 if(fXSType == fEmIncreasing) { 997 if(fXSType == fEmIncreasing) {
689 if(e*invLambdaFactor < mfpKinEnergy) { << 998 if(e/lambdaFactor < mfpKinEnergy) {
690 preStepLambda = GetLambdaForScaledEnergy << 999 mfpKinEnergy = e;
691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 1000 preStepLambda = GetLambdaForScaledEnergy(e, loge);
692 } <<
693 <<
694 // cross section has one peak <<
695 } else if(fXSType == fEmOnePeak) { <<
696 const G4double epeak = (*theEnergyOfCrossS <<
697 if(e <= epeak) { <<
698 if(e*invLambdaFactor < mfpKinEnergy) { <<
699 preStepLambda = GetLambdaForScaledEner <<
700 mfpKinEnergy = (preStepLambda > 0.0) ? <<
701 } <<
702 } else if(e < mfpKinEnergy) { <<
703 const G4double e1 = std::max(epeak, e*la <<
704 mfpKinEnergy = e1; <<
705 preStepLambda = GetLambdaForScaledEnergy <<
706 } 1001 }
707 1002
708 // cross section has more than one peaks << 1003 // cross section has two peaks
709 } else if(fXSType == fEmTwoPeaks) { 1004 } else if(fXSType == fEmTwoPeaks) {
710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple 1005 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCoupleIndex];
711 const G4double e1peak = xs->e1peak; 1006 const G4double e1peak = xs->e1peak;
712 1007
713 // below the 1st peak 1008 // below the 1st peak
714 if(e <= e1peak) { 1009 if(e <= e1peak) {
715 if(e*invLambdaFactor < mfpKinEnergy) { << 1010 if(e/lambdaFactor < mfpKinEnergy) {
716 preStepLambda = GetLambdaForScaledEner << 1011 mfpKinEnergy = e;
717 mfpKinEnergy = (preStepLambda > 0.0) ? << 1012 preStepLambda = GetLambdaForScaledEnergy(e, loge);
718 } 1013 }
719 return; 1014 return;
720 } 1015 }
721 const G4double e1deep = xs->e1deep; 1016 const G4double e1deep = xs->e1deep;
722 // above the 1st peak, below the deep 1017 // above the 1st peak, below the deep
723 if(e <= e1deep) { 1018 if(e <= e1deep) {
724 if(mfpKinEnergy >= e1deep || e <= mfpKin 1019 if(mfpKinEnergy >= e1deep || e <= mfpKinEnergy) {
725 const G4double e1 = std::max(e1peak, e 1020 const G4double e1 = std::max(e1peak, e*lambdaFactor);
726 mfpKinEnergy = e1; <<
727 preStepLambda = GetLambdaForScaledEner 1021 preStepLambda = GetLambdaForScaledEnergy(e1);
>> 1022 mfpKinEnergy = e1;
728 } 1023 }
729 return; 1024 return;
730 } 1025 }
731 const G4double e2peak = xs->e2peak; 1026 const G4double e2peak = xs->e2peak;
732 // above the deep, below 2nd peak 1027 // above the deep, below 2nd peak
733 if(e <= e2peak) { 1028 if(e <= e2peak) {
734 if(e*invLambdaFactor < mfpKinEnergy) { << 1029 if(e/lambdaFactor < mfpKinEnergy) {
735 mfpKinEnergy = e; 1030 mfpKinEnergy = e;
736 preStepLambda = GetLambdaForScaledEner << 1031 preStepLambda = GetLambdaForScaledEnergy(e, loge);
737 } 1032 }
738 return; 1033 return;
739 } 1034 }
740 const G4double e2deep = xs->e2deep; 1035 const G4double e2deep = xs->e2deep;
741 // above the 2nd peak, below the deep 1036 // above the 2nd peak, below the deep
742 if(e <= e2deep) { 1037 if(e <= e2deep) {
743 if(mfpKinEnergy >= e2deep || e <= mfpKin 1038 if(mfpKinEnergy >= e2deep || e <= mfpKinEnergy) {
744 const G4double e1 = std::max(e2peak, e 1039 const G4double e1 = std::max(e2peak, e*lambdaFactor);
745 mfpKinEnergy = e1; <<
746 preStepLambda = GetLambdaForScaledEner 1040 preStepLambda = GetLambdaForScaledEnergy(e1);
>> 1041 mfpKinEnergy = e1;
747 } 1042 }
748 return; 1043 return;
749 } 1044 }
750 const G4double e3peak = xs->e3peak; <<
751 // above the deep, below 3d peak 1045 // above the deep, below 3d peak
752 if(e <= e3peak) { << 1046 if(e/lambdaFactor < mfpKinEnergy) {
753 if(e*invLambdaFactor < mfpKinEnergy) { << 1047 mfpKinEnergy = e;
754 mfpKinEnergy = e; << 1048 preStepLambda = GetLambdaForScaledEnergy(e, loge);
755 preStepLambda = GetLambdaForScaledEner <<
756 } <<
757 return; <<
758 } <<
759 // above 3d peak <<
760 if(e <= mfpKinEnergy) { <<
761 const G4double e1 = std::max(e3peak, e*l <<
762 mfpKinEnergy = e1; <<
763 preStepLambda = GetLambdaForScaledEnergy <<
764 } 1049 }
>> 1050
765 // integral method is not used 1051 // integral method is not used
766 } else { 1052 } else {
767 preStepLambda = GetLambdaForScaledEnergy(e << 1053 preStepLambda = GetLambdaForScaledEnergy(e, loge);
768 } 1054 }
769 } 1055 }
770 1056
771 //....oooOO0OOooo........oooOO0OOooo........oo 1057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
772 1058
773 G4VParticleChange* G4VEnergyLossProcess::Along 1059 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
774 1060 const G4Step& step)
775 { 1061 {
776 fParticleChange.InitializeForAlongStep(track 1062 fParticleChange.InitializeForAlongStep(track);
777 // The process has range table - calculate e 1063 // The process has range table - calculate energy loss
778 if(!isIonisation || !currentModel->IsActive( 1064 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) {
779 return &fParticleChange; 1065 return &fParticleChange;
780 } 1066 }
781 1067
>> 1068 // Get the actual (true) Step length
782 G4double length = step.GetStepLength(); 1069 G4double length = step.GetStepLength();
>> 1070 if(length <= 0.0) { return &fParticleChange; }
783 G4double eloss = 0.0; 1071 G4double eloss = 0.0;
784 1072
785 /* << 1073 /*
786 if(-1 < verboseLevel) { 1074 if(-1 < verboseLevel) {
787 const G4ParticleDefinition* d = track.GetP 1075 const G4ParticleDefinition* d = track.GetParticleDefinition();
788 G4cout << "AlongStepDoIt for " 1076 G4cout << "AlongStepDoIt for "
789 << GetProcessName() << " and partic << 1077 << GetProcessName() << " and particle "
790 << " eScaled(MeV)=" << preStepScal << 1078 << d->GetParticleName()
791 << " range(mm)=" << fRange/mm << " << 1079 << " eScaled(MeV)= " << preStepScaledEnergy/MeV
792 << " rf=" << reduceFactor << " q^ << 1080 << " range(mm)= " << fRange/mm
793 << " md=" << d->GetPDGMass() << " << 1081 << " s(mm)= " << length/mm
794 << " " << track.GetMaterial()->Get << 1082 << " rf= " << reduceFactor
>> 1083 << " q^2= " << chargeSqRatio
>> 1084 << " md= " << d->GetPDGMass()
>> 1085 << " status= " << track.GetTrackStatus()
>> 1086 << " " << track.GetMaterial()->GetName()
>> 1087 << G4endl;
795 } 1088 }
796 */ 1089 */
>> 1090
797 const G4DynamicParticle* dynParticle = track 1091 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
798 1092
799 // define new weight for primary and seconda 1093 // define new weight for primary and secondaries
800 G4double weight = fParticleChange.GetParentW 1094 G4double weight = fParticleChange.GetParentWeight();
801 if(weightFlag) { 1095 if(weightFlag) {
802 weight /= biasFactor; 1096 weight /= biasFactor;
803 fParticleChange.ProposeWeight(weight); 1097 fParticleChange.ProposeWeight(weight);
804 } 1098 }
805 1099
806 // stopping, check actual range and kinetic << 1100 // stopping
807 if (length >= fRange || preStepKinEnergy <= 1101 if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) {
808 eloss = preStepKinEnergy; 1102 eloss = preStepKinEnergy;
809 if (useDeexcitation) { 1103 if (useDeexcitation) {
810 atomDeexcitation->AlongStepDeexcitation( 1104 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
811 << 1105 eloss, currentCoupleIndex);
812 if(scTracks.size() > 0) { FillSecondarie 1106 if(scTracks.size() > 0) { FillSecondariesAlongStep(weight); }
813 eloss = std::max(eloss, 0.0); 1107 eloss = std::max(eloss, 0.0);
814 } 1108 }
815 fParticleChange.SetProposedKineticEnergy(0 1109 fParticleChange.SetProposedKineticEnergy(0.0);
816 fParticleChange.ProposeLocalEnergyDeposit( 1110 fParticleChange.ProposeLocalEnergyDeposit(eloss);
817 return &fParticleChange; 1111 return &fParticleChange;
818 } 1112 }
819 // zero step length with non-zero range << 1113 //G4cout << theDEDXTable << " idx= " << basedCoupleIndex
820 if(length <= 0.0) { return &fParticleChange; << 1114 // << " " << GetProcessName() << " "<< currentMaterial->GetName()<<G4endl;
821 << 1115 //if(particle->GetParticleName() == "e-")G4cout << (*theDEDXTable) <<G4endl;
822 // Short step 1116 // Short step
823 eloss = length*GetDEDXForScaledEnergy(preSte 1117 eloss = length*GetDEDXForScaledEnergy(preStepScaledEnergy,
824 LogSca << 1118 preStepLogScaledEnergy);
825 /* << 1119
826 G4cout << "##### Short STEP: eloss= " << elo << 1120 //G4cout << "Short STEP: eloss= " << eloss << G4endl;
827 << " Escaled=" << preStepScaledEnergy << 1121
828 << " R=" << fRange <<
829 << " L=" << length <<
830 << " fFactor=" << fFactor << " minE=" << mi <<
831 << " idxBase=" << basedCoupleIndex << G4end <<
832 */ <<
833 // Long step 1122 // Long step
834 if(eloss > preStepKinEnergy*linLossLimit) { 1123 if(eloss > preStepKinEnergy*linLossLimit) {
835 1124
836 const G4double x = (fRange - length)/reduc << 1125 G4double x = (fRange - length)/reduceFactor;
837 const G4double de = preStepKinEnergy - Sca << 1126 //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl;
838 if(de > 0.0) { eloss = de; } << 1127 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
>> 1128
839 /* 1129 /*
840 if(-1 < verboseLevel) 1130 if(-1 < verboseLevel)
841 G4cout << " Long STEP: rPre(mm)=" << 1131 G4cout << "Long STEP: rPre(mm)= "
842 << GetScaledRangeForScaledEnergy( 1132 << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
843 << " x(mm)=" << x/mm << 1133 << " rPost(mm)= " << x/mm
844 << " eloss(MeV)=" << eloss/MeV << 1134 << " ePre(MeV)= " << preStepScaledEnergy/MeV
845 << " rFactor=" << reduceFactor << 1135 << " eloss(MeV)= " << eloss/MeV
846 << " massRatio=" << massRatio << 1136 << " eloss0(MeV)= "
>> 1137 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
>> 1138 << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
847 << G4endl; 1139 << G4endl;
848 */ 1140 */
849 } 1141 }
850 1142
851 /* 1143 /*
>> 1144 G4double eloss0 = eloss;
852 if(-1 < verboseLevel ) { 1145 if(-1 < verboseLevel ) {
853 G4cout << "Before fluct: eloss(MeV)= " << 1146 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
854 << " e-eloss= " << preStepKinEnergy 1147 << " e-eloss= " << preStepKinEnergy-eloss
855 << " step(mm)= " << length/mm << " << 1148 << " step(mm)= " << length/mm
856 << " fluct= " << lossFluctuationFla << 1149 << " range(mm)= " << fRange/mm
>> 1150 << " fluct= " << lossFluctuationFlag
>> 1151 << G4endl;
857 } 1152 }
858 */ 1153 */
859 1154
860 const G4double cut = (*theCuts)[currentCoupl 1155 const G4double cut = (*theCuts)[currentCoupleIndex];
861 G4double esec = 0.0; 1156 G4double esec = 0.0;
862 1157
863 // Corrections, which cannot be tabulated 1158 // Corrections, which cannot be tabulated
864 if(isIon) { 1159 if(isIon) {
865 currentModel->CorrectionsAlongStep(current 1160 currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
866 length, 1161 length, eloss);
867 eloss = std::max(eloss, 0.0); 1162 eloss = std::max(eloss, 0.0);
868 } 1163 }
869 1164
870 // Sample fluctuations if not full energy lo 1165 // Sample fluctuations if not full energy loss
871 if(eloss >= preStepKinEnergy) { 1166 if(eloss >= preStepKinEnergy) {
872 eloss = preStepKinEnergy; 1167 eloss = preStepKinEnergy;
873 1168
874 } else if (lossFluctuationFlag) { 1169 } else if (lossFluctuationFlag) {
875 const G4double tmax = currentModel->MaxSec 1170 const G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
876 const G4double tcut = std::min(cut, tmax); 1171 const G4double tcut = std::min(cut, tmax);
877 G4VEmFluctuationModel* fluc = currentModel 1172 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
878 eloss = fluc->SampleFluctuations(currentCo 1173 eloss = fluc->SampleFluctuations(currentCouple,dynParticle,
879 tcut, tma 1174 tcut, tmax, length, eloss);
880 /* 1175 /*
881 if(-1 < verboseLevel) 1176 if(-1 < verboseLevel)
882 G4cout << "After fluct: eloss(MeV)= " << 1177 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
883 << " fluc= " << (eloss-eloss0)/Me 1178 << " fluc= " << (eloss-eloss0)/MeV
884 << " ChargeSqRatio= " << chargeSq 1179 << " ChargeSqRatio= " << chargeSqRatio
885 << " massRatio= " << massRatio << << 1180 << " massRatio= " << massRatio
>> 1181 << " tmax= " << tmax
>> 1182 << G4endl;
886 */ 1183 */
887 } 1184 }
888 1185
889 // deexcitation 1186 // deexcitation
890 if (useDeexcitation) { 1187 if (useDeexcitation) {
891 G4double esecfluo = preStepKinEnergy; 1188 G4double esecfluo = preStepKinEnergy;
892 G4double de = esecfluo; 1189 G4double de = esecfluo;
>> 1190 //G4double eloss0 = eloss;
>> 1191 /*
>> 1192 G4cout << "### 1: E(keV)= " << preStepKinEnergy/keV
>> 1193 << " Efluomax(keV)= " << de/keV
>> 1194 << " Eloss(keV)= " << eloss/keV << G4endl;
>> 1195 */
893 atomDeexcitation->AlongStepDeexcitation(sc 1196 atomDeexcitation->AlongStepDeexcitation(scTracks, step,
894 de << 1197 de, currentCoupleIndex);
895 1198
896 // sum of de-excitation energies 1199 // sum of de-excitation energies
897 esecfluo -= de; 1200 esecfluo -= de;
898 1201
899 // subtracted from energy loss 1202 // subtracted from energy loss
900 if(eloss >= esecfluo) { 1203 if(eloss >= esecfluo) {
901 esec += esecfluo; 1204 esec += esecfluo;
902 eloss -= esecfluo; 1205 eloss -= esecfluo;
903 } else { 1206 } else {
904 esec += esecfluo; 1207 esec += esecfluo;
905 eloss = 0.0; 1208 eloss = 0.0;
906 } 1209 }
>> 1210 /*
>> 1211 if(esecfluo > 0.0) {
>> 1212 G4cout << "### 2: E(keV)= " << preStepKinEnergy/keV
>> 1213 << " Esec(keV)= " << esec/keV
>> 1214 << " Esecf(kV)= " << esecfluo/keV
>> 1215 << " Eloss0(kV)= " << eloss0/keV
>> 1216 << " Eloss(keV)= " << eloss/keV
>> 1217 << G4endl;
>> 1218 }
>> 1219 */
907 } 1220 }
908 if(nullptr != subcutProducer && IsRegionForC 1221 if(nullptr != subcutProducer && IsRegionForCubcutProcessor(track)) {
909 subcutProducer->SampleSecondaries(step, sc 1222 subcutProducer->SampleSecondaries(step, scTracks, eloss, cut);
910 } 1223 }
911 // secondaries from atomic de-excitation and 1224 // secondaries from atomic de-excitation and subcut
912 if(!scTracks.empty()) { FillSecondariesAlong 1225 if(!scTracks.empty()) { FillSecondariesAlongStep(weight); }
913 1226
914 // Energy balance 1227 // Energy balance
915 G4double finalT = preStepKinEnergy - eloss - 1228 G4double finalT = preStepKinEnergy - eloss - esec;
916 if (finalT <= lowestKinEnergy) { 1229 if (finalT <= lowestKinEnergy) {
917 eloss += finalT; 1230 eloss += finalT;
918 finalT = 0.0; 1231 finalT = 0.0;
919 } else if(isIon) { 1232 } else if(isIon) {
920 fParticleChange.SetProposedCharge( 1233 fParticleChange.SetProposedCharge(
921 currentModel->GetParticleCharge(track.Ge 1234 currentModel->GetParticleCharge(track.GetParticleDefinition(),
922 currentM 1235 currentMaterial,finalT));
923 } 1236 }
924 eloss = std::max(eloss, 0.0); 1237 eloss = std::max(eloss, 0.0);
925 1238
926 fParticleChange.SetProposedKineticEnergy(fin 1239 fParticleChange.SetProposedKineticEnergy(finalT);
927 fParticleChange.ProposeLocalEnergyDeposit(el 1240 fParticleChange.ProposeLocalEnergyDeposit(eloss);
928 /* 1241 /*
929 if(-1 < verboseLevel) { 1242 if(-1 < verboseLevel) {
930 G4double del = finalT + eloss + esec - pre 1243 G4double del = finalT + eloss + esec - preStepKinEnergy;
931 G4cout << "Final value eloss(MeV)= " << el 1244 G4cout << "Final value eloss(MeV)= " << eloss/MeV
932 << " preStepKinEnergy= " << preStep 1245 << " preStepKinEnergy= " << preStepKinEnergy
933 << " postStepKinEnergy= " << finalT 1246 << " postStepKinEnergy= " << finalT
934 << " de(keV)= " << del/keV 1247 << " de(keV)= " << del/keV
935 << " lossFlag= " << lossFluctuation 1248 << " lossFlag= " << lossFluctuationFlag
936 << " status= " << track.GetTrackSt 1249 << " status= " << track.GetTrackStatus()
937 << G4endl; 1250 << G4endl;
938 } 1251 }
939 */ 1252 */
940 return &fParticleChange; 1253 return &fParticleChange;
941 } 1254 }
942 1255
943 //....oooOO0OOooo........oooOO0OOooo........oo 1256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
944 1257
945 void G4VEnergyLossProcess::FillSecondariesAlon 1258 void G4VEnergyLossProcess::FillSecondariesAlongStep(G4double wt)
946 { 1259 {
947 const std::size_t n0 = scTracks.size(); << 1260 const G4int n0 = scTracks.size();
948 G4double weight = wt; 1261 G4double weight = wt;
949 // weight may be changed by biasing manager 1262 // weight may be changed by biasing manager
950 if(biasManager) { 1263 if(biasManager) {
951 if(biasManager->SecondaryBiasingRegion((G4 << 1264 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
952 weight *= 1265 weight *=
953 biasManager->ApplySecondaryBiasing(scT << 1266 biasManager->ApplySecondaryBiasing(scTracks, currentCoupleIndex);
954 } 1267 }
955 } 1268 }
956 1269
957 // fill secondaries 1270 // fill secondaries
958 const std::size_t n = scTracks.size(); << 1271 const G4int n = scTracks.size();
959 fParticleChange.SetNumberOfSecondaries((G4in << 1272 fParticleChange.SetNumberOfSecondaries(n);
960 1273
961 for(std::size_t i=0; i<n; ++i) { << 1274 for(G4int i=0; i<n; ++i) {
962 G4Track* t = scTracks[i]; 1275 G4Track* t = scTracks[i];
963 if(nullptr != t) { 1276 if(nullptr != t) {
964 t->SetWeight(weight); 1277 t->SetWeight(weight);
965 pParticleChange->AddSecondary(t); 1278 pParticleChange->AddSecondary(t);
966 G4int pdg = t->GetDefinition()->GetPDGEn << 1279 if(i >= n0) { t->SetCreatorModelID(biasID); }
967 if (i < n0) { << 1280 //G4cout << "Secondary(along step) has weight " << t->GetWeight()
968 if (pdg == 22) { << 1281 //<< ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
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 } 1282 }
979 } 1283 }
980 scTracks.clear(); 1284 scTracks.clear();
981 } 1285 }
982 1286
983 //....oooOO0OOooo........oooOO0OOooo........oo 1287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
984 1288
985 G4VParticleChange* G4VEnergyLossProcess::PostS 1289 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
986 1290 const G4Step& step)
987 { 1291 {
988 // clear number of interaction lengths in an << 1292 // In all cases clear number of interaction lengths
989 theNumberOfInteractionLengthLeft = -1.0; 1293 theNumberOfInteractionLengthLeft = -1.0;
990 mfpKinEnergy = DBL_MAX; 1294 mfpKinEnergy = DBL_MAX;
991 1295
992 fParticleChange.InitializeForPostStep(track) 1296 fParticleChange.InitializeForPostStep(track);
993 const G4double finalT = track.GetKineticEner 1297 const G4double finalT = track.GetKineticEnergy();
994 1298
995 const G4double postStepScaledEnergy = finalT 1299 const G4double postStepScaledEnergy = finalT*massRatio;
996 SelectModel(postStepScaledEnergy); 1300 SelectModel(postStepScaledEnergy);
997 1301
998 if(!currentModel->IsActive(postStepScaledEne 1302 if(!currentModel->IsActive(postStepScaledEnergy)) {
999 return &fParticleChange; 1303 return &fParticleChange;
1000 } 1304 }
1001 /* 1305 /*
1002 if(1 < verboseLevel) { << 1306 if(-1 < verboseLevel) {
1003 G4cout<<GetProcessName()<<" PostStepDoIt: << 1307 G4cout << GetProcessName()
>> 1308 << "::PostStepDoIt: E(MeV)= " << finalT/MeV
>> 1309 << G4endl;
1004 } 1310 }
1005 */ 1311 */
>> 1312
1006 // forced process - should happen only once 1313 // forced process - should happen only once per track
1007 if(biasFlag) { 1314 if(biasFlag) {
1008 if(biasManager->ForcedInteractionRegion(( << 1315 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
1009 biasFlag = false; 1316 biasFlag = false;
1010 } 1317 }
1011 } 1318 }
>> 1319
1012 const G4DynamicParticle* dp = track.GetDyna 1320 const G4DynamicParticle* dp = track.GetDynamicParticle();
1013 1321
1014 // Integral approach 1322 // Integral approach
1015 if (fXSType != fEmNoIntegral) { 1323 if (fXSType != fEmNoIntegral) {
1016 const G4double logFinalT = dp->GetLogKine 1324 const G4double logFinalT = dp->GetLogKineticEnergy();
1017 G4double lx = GetLambdaForScaledEnergy(po 1325 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy,
1018 lo 1326 logFinalT + logMassRatio);
1019 lx = std::max(lx, 0.0); 1327 lx = std::max(lx, 0.0);
1020 1328
>> 1329 // cache cross section useful for the false interaction
>> 1330 const G4double lg = preStepLambda;
>> 1331 if(postStepScaledEnergy < mfpKinEnergy) {
>> 1332 mfpKinEnergy = postStepScaledEnergy;
>> 1333 preStepLambda = lx;
>> 1334 }
>> 1335 /*
>> 1336 if(preStepLambda<lx && 1 < verboseLevel) {
>> 1337 G4cout << "WARNING: for " << particle->GetParticleName()
>> 1338 << " and " << GetProcessName()
>> 1339 << " E(MeV)= " << finalT/MeV
>> 1340 << " preLambda= " << preStepLambda
>> 1341 << " < " << lx << " (postLambda) "
>> 1342 << G4endl;
>> 1343 }
>> 1344 */
1021 // if both lg and lx are zero then no int 1345 // if both lg and lx are zero then no interaction
1022 if(preStepLambda*G4UniformRand() >= lx) { << 1346 if(lg*G4UniformRand() >= lx) {
1023 return &fParticleChange; 1347 return &fParticleChange;
1024 } 1348 }
1025 } 1349 }
1026 <<
1027 // define new weight for primary and second 1350 // define new weight for primary and secondaries
1028 G4double weight = fParticleChange.GetParent 1351 G4double weight = fParticleChange.GetParentWeight();
1029 if(weightFlag) { 1352 if(weightFlag) {
1030 weight /= biasFactor; 1353 weight /= biasFactor;
1031 fParticleChange.ProposeWeight(weight); 1354 fParticleChange.ProposeWeight(weight);
1032 } 1355 }
1033 1356
1034 const G4double tcut = (*theCuts)[currentCou 1357 const G4double tcut = (*theCuts)[currentCoupleIndex];
1035 1358
1036 // sample secondaries 1359 // sample secondaries
1037 secParticles.clear(); 1360 secParticles.clear();
>> 1361 //G4cout<< "@@@ Eprimary= "<<dynParticle->GetKineticEnergy()/MeV
>> 1362 // << " cut= " << tcut/MeV << G4endl;
1038 currentModel->SampleSecondaries(&secParticl 1363 currentModel->SampleSecondaries(&secParticles, currentCouple, dp, tcut);
1039 1364
1040 const G4int num0 = (G4int)secParticles.size << 1365 const G4int num0 = secParticles.size();
1041 1366
1042 // bremsstrahlung splitting or Russian roul 1367 // bremsstrahlung splitting or Russian roulette
1043 if(biasManager) { 1368 if(biasManager) {
1044 if(biasManager->SecondaryBiasingRegion((G << 1369 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
1045 G4double eloss = 0.0; 1370 G4double eloss = 0.0;
1046 weight *= biasManager->ApplySecondaryBi 1371 weight *= biasManager->ApplySecondaryBiasing(
1047 secPart 1372 secParticles,
1048 track, 1373 track, currentModel,
1049 &fParti 1374 &fParticleChange, eloss,
1050 (G4int) << 1375 currentCoupleIndex, tcut,
1051 step.Ge 1376 step.GetPostStepPoint()->GetSafety());
1052 if(eloss > 0.0) { 1377 if(eloss > 0.0) {
1053 eloss += fParticleChange.GetLocalEner 1378 eloss += fParticleChange.GetLocalEnergyDeposit();
1054 fParticleChange.ProposeLocalEnergyDep 1379 fParticleChange.ProposeLocalEnergyDeposit(eloss);
1055 } 1380 }
1056 } 1381 }
1057 } 1382 }
1058 1383
1059 // save secondaries 1384 // save secondaries
1060 const G4int num = (G4int)secParticles.size( << 1385 const G4int num = secParticles.size();
1061 if(num > 0) { 1386 if(num > 0) {
1062 1387
1063 fParticleChange.SetNumberOfSecondaries(nu 1388 fParticleChange.SetNumberOfSecondaries(num);
1064 G4double time = track.GetGlobalTime(); 1389 G4double time = track.GetGlobalTime();
1065 1390
1066 G4int n1(0), n2(0); 1391 G4int n1(0), n2(0);
1067 if(num0 > mainSecondaries) { 1392 if(num0 > mainSecondaries) {
1068 currentModel->FillNumberOfSecondaries(n 1393 currentModel->FillNumberOfSecondaries(n1, n2);
1069 } 1394 }
1070 1395
1071 for (G4int i=0; i<num; ++i) { 1396 for (G4int i=0; i<num; ++i) {
1072 if(nullptr != secParticles[i]) { 1397 if(nullptr != secParticles[i]) {
1073 G4Track* t = new G4Track(secParticles 1398 G4Track* t = new G4Track(secParticles[i], time, track.GetPosition());
1074 t->SetTouchableHandle(track.GetToucha 1399 t->SetTouchableHandle(track.GetTouchableHandle());
1075 if (biasManager) { 1400 if (biasManager) {
1076 t->SetWeight(weight * biasManager-> 1401 t->SetWeight(weight * biasManager->GetWeight(i));
1077 } else { 1402 } else {
1078 t->SetWeight(weight); 1403 t->SetWeight(weight);
1079 } 1404 }
1080 if(i < num0) { 1405 if(i < num0) {
1081 t->SetCreatorModelID(secID); 1406 t->SetCreatorModelID(secID);
1082 } else if(i < num0 + n1) { 1407 } else if(i < num0 + n1) {
1083 t->SetCreatorModelID(tripletID); 1408 t->SetCreatorModelID(tripletID);
1084 } else { 1409 } else {
1085 t->SetCreatorModelID(biasID); 1410 t->SetCreatorModelID(biasID);
1086 } 1411 }
1087 1412
1088 //G4cout << "Secondary(post step) has 1413 //G4cout << "Secondary(post step) has weight " << t->GetWeight()
1089 // << ", kenergy " << t->GetKin 1414 // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV"
1090 // << " time= " << time/ns << " 1415 // << " time= " << time/ns << " ns " << G4endl;
1091 pParticleChange->AddSecondary(t); 1416 pParticleChange->AddSecondary(t);
1092 } 1417 }
1093 } 1418 }
1094 } 1419 }
1095 1420
1096 if(0.0 == fParticleChange.GetProposedKineti 1421 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
1097 fAlive == fParticleChange.GetTrackStatus 1422 fAlive == fParticleChange.GetTrackStatus()) {
1098 if(particle->GetProcessManager()->GetAtRe 1423 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
1099 { fParticleChange.ProposeTrackStatus 1424 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
1100 else { fParticleChange.ProposeTrackStatus 1425 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
1101 } 1426 }
1102 1427
1103 /* 1428 /*
1104 if(-1 < verboseLevel) { 1429 if(-1 < verboseLevel) {
1105 G4cout << "::PostStepDoIt: Sample seconda 1430 G4cout << "::PostStepDoIt: Sample secondary; Efin= "
1106 << fParticleChange.GetProposedKineticEner 1431 << fParticleChange.GetProposedKineticEnergy()/MeV
1107 << " MeV; model= (" << currentMode 1432 << " MeV; model= (" << currentModel->LowEnergyLimit()
1108 << ", " << currentModel->HighEner 1433 << ", " << currentModel->HighEnergyLimit() << ")"
1109 << " preStepLambda= " << preStepL 1434 << " preStepLambda= " << preStepLambda
1110 << " dir= " << track.GetMomentumD 1435 << " dir= " << track.GetMomentumDirection()
1111 << " status= " << track.GetTrackS 1436 << " status= " << track.GetTrackStatus()
1112 << G4endl; 1437 << G4endl;
1113 } 1438 }
1114 */ 1439 */
1115 return &fParticleChange; 1440 return &fParticleChange;
1116 } 1441 }
1117 1442
1118 //....oooOO0OOooo........oooOO0OOooo........o 1443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1119 1444
1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl 1445 G4bool G4VEnergyLossProcess::StorePhysicsTable(
1121 const G4ParticleDefinition* part, cons << 1446 const G4ParticleDefinition* part, const G4String& directory,
>> 1447 G4bool ascii)
1122 { 1448 {
1123 if (!isMaster || nullptr != baseParticle || << 1449 G4bool res = true;
1124 for(std::size_t i=0; i<7; ++i) { << 1450 //G4cout << "G4VEnergyLossProcess::StorePhysicsTable: " << part->GetParticleName()
1125 // ionisation table only for ionisation p << 1451 // << " " << directory << " " << ascii << G4endl;
1126 if (nullptr == theData->Table(i) || (!isI << 1452 if (!isMaster || baseParticle || part != particle ) return res;
1127 continue; << 1453
1128 } << 1454 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
1129 if (-1 < verboseLevel) { << 1455 {res = false;}
1130 G4cout << "G4VEnergyLossProcess::StoreP << 1456
1131 << " " << particle->GetParticleName() << 1457 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
1132 << " " << GetProcessName() << 1458 {res = false;}
1133 << " " << tnames[i] << " " << theDat << 1459
1134 } << 1460 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
1135 if (!G4EmTableUtil::StoreTable(this, part << 1461 {res = false;}
1136 dir, tnames[i], verboseLevel, asci << 1462
1137 return false; << 1463 if(isIonisation &&
1138 } << 1464 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
1139 } << 1465 {res = false;}
1140 return true; << 1466
>> 1467 if(isIonisation &&
>> 1468 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
>> 1469 {res = false;}
>> 1470
>> 1471 if(isIonisation &&
>> 1472 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
>> 1473 {res = false;}
>> 1474
>> 1475 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
>> 1476 {res = false;}
>> 1477
>> 1478 return res;
1141 } 1479 }
1142 1480
1143 //....oooOO0OOooo........oooOO0OOooo........o 1481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1144 1482
1145 G4bool 1483 G4bool
1146 G4VEnergyLossProcess::RetrievePhysicsTable(co 1484 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
1147 co << 1485 const G4String& directory,
>> 1486 G4bool ascii)
1148 { 1487 {
1149 if (!isMaster || nullptr != baseParticle || << 1488 G4bool res = true;
1150 for(std::size_t i=0; i<7; ++i) { << 1489 if (!isMaster) return res;
1151 // ionisation table only for ionisation p << 1490 const G4String& particleName = part->GetParticleName();
1152 if (!isIonisation && 1 == i) { continue; << 1491
1153 if(!G4EmTableUtil::RetrieveTable(this, pa << 1492 if(1 < verboseLevel) {
1154 verboseL << 1493 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
1155 return false; << 1494 << particleName << " and process " << GetProcessName()
>> 1495 << "; tables_are_built= " << tablesAreBuilt
>> 1496 << G4endl;
>> 1497 }
>> 1498 if(particle == part) {
>> 1499
>> 1500 if(nullptr == baseParticle) {
>> 1501
>> 1502 G4bool fpi = true;
>> 1503 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
>> 1504 { fpi = false; }
>> 1505
>> 1506 // ionisation table keeps individual dEdx and not sum of sub-processes
>> 1507 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false))
>> 1508 { fpi = false; }
>> 1509
>> 1510 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
>> 1511 { res = false; }
>> 1512
>> 1513 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,
>> 1514 "DEDXnr",false))
>> 1515 { res = false; }
>> 1516
>> 1517 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,
>> 1518 "CSDARange",false))
>> 1519 { res = false; }
>> 1520
>> 1521 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,
>> 1522 "InverseRange",fpi))
>> 1523 { res = false; }
>> 1524
>> 1525 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
>> 1526 { res = false; }
1156 } 1527 }
1157 } 1528 }
>> 1529 return res;
>> 1530 }
>> 1531
>> 1532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
>> 1533
>> 1534 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
>> 1535 G4PhysicsTable* aTable, G4bool ascii,
>> 1536 const G4String& directory,
>> 1537 const G4String& tname)
>> 1538 {
>> 1539 G4bool res = true;
>> 1540 if (nullptr != aTable) {
>> 1541 const G4String& name = GetPhysicsTableFileName(part, directory, tname, ascii);
>> 1542 if ( aTable->StorePhysicsTable(name,ascii) ) {
>> 1543 if (0 < verboseLevel) G4cout << "Stored: " << name << G4endl;
>> 1544 } else {
>> 1545 res = false;
>> 1546 G4cout << "Fail to store: " << name << G4endl;
>> 1547 }
>> 1548 }
>> 1549 return res;
>> 1550 }
>> 1551
>> 1552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
>> 1553
>> 1554 G4bool
>> 1555 G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
>> 1556 G4PhysicsTable* aTable,
>> 1557 G4bool ascii,
>> 1558 const G4String& directory,
>> 1559 const G4String& tname,
>> 1560 G4bool mandatory)
>> 1561 {
>> 1562 G4bool isRetrieved = false;
>> 1563 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
>> 1564 if(nullptr != aTable) {
>> 1565 if(aTable->ExistPhysicsTable(filename)) {
>> 1566 if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii,spline)) {
>> 1567 isRetrieved = true;
>> 1568 if(spline) {
>> 1569 for(auto & v : *aTable) {
>> 1570 if(nullptr != v) { v->FillSecondDerivatives(); }
>> 1571 }
>> 1572 }
>> 1573 if (0 < verboseLevel) {
>> 1574 G4cout << tname << " table for " << part->GetParticleName()
>> 1575 << " is Retrieved from <" << filename << ">"
>> 1576 << G4endl;
>> 1577 }
>> 1578 }
>> 1579 }
>> 1580 }
>> 1581 if(mandatory && !isRetrieved) {
>> 1582 if(0 < verboseLevel) {
>> 1583 G4cout << tname << " table for " << part->GetParticleName()
>> 1584 << " from file <"
>> 1585 << filename << "> is not Retrieved"
>> 1586 << G4endl;
>> 1587 }
>> 1588 return false;
>> 1589 }
1158 return true; 1590 return true;
1159 } 1591 }
1160 1592
1161 //....oooOO0OOooo........oooOO0OOooo........o 1593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1162 1594
1163 G4double G4VEnergyLossProcess::GetDEDXDispers 1595 G4double G4VEnergyLossProcess::GetDEDXDispersion(
1164 const G4Mat 1596 const G4MaterialCutsCouple *couple,
1165 const G4Dyn 1597 const G4DynamicParticle* dp,
1166 G4dou 1598 G4double length)
1167 { 1599 {
1168 DefineMaterial(couple); 1600 DefineMaterial(couple);
1169 G4double ekin = dp->GetKineticEnergy(); 1601 G4double ekin = dp->GetKineticEnergy();
1170 SelectModel(ekin*massRatio); 1602 SelectModel(ekin*massRatio);
1171 G4double tmax = currentModel->MaxSecondaryK 1603 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
1172 G4double tcut = std::min(tmax,(*theCuts)[cu 1604 G4double tcut = std::min(tmax,(*theCuts)[currentCoupleIndex]);
1173 G4double d = 0.0; 1605 G4double d = 0.0;
1174 G4VEmFluctuationModel* fm = currentModel->G 1606 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
1175 if(nullptr != fm) { d = fm->Dispersion(curr 1607 if(nullptr != fm) { d = fm->Dispersion(currentMaterial,dp,tcut,tmax,length); }
1176 return d; 1608 return d;
1177 } 1609 }
1178 1610
1179 //....oooOO0OOooo........oooOO0OOooo........o 1611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1180 1612
1181 G4double 1613 G4double
1182 G4VEnergyLossProcess::CrossSectionPerVolume(G 1614 G4VEnergyLossProcess::CrossSectionPerVolume(G4double kineticEnergy,
1183 c 1615 const G4MaterialCutsCouple* couple,
1184 G 1616 G4double logKineticEnergy)
1185 { 1617 {
1186 // Cross section per volume is calculated 1618 // Cross section per volume is calculated
1187 DefineMaterial(couple); 1619 DefineMaterial(couple);
1188 G4double cross = 0.0; 1620 G4double cross = 0.0;
1189 if (nullptr != theLambdaTable) { 1621 if (nullptr != theLambdaTable) {
1190 cross = GetLambdaForScaledEnergy(kineticE 1622 cross = GetLambdaForScaledEnergy(kineticEnergy * massRatio,
1191 logKinet 1623 logKineticEnergy + logMassRatio);
1192 } else { 1624 } else {
1193 SelectModel(kineticEnergy*massRatio); 1625 SelectModel(kineticEnergy*massRatio);
1194 cross = (!baseMat) ? biasFactor : biasFac 1626 cross = (!baseMat) ? biasFactor : biasFactor*(*theDensityFactor)[currentCoupleIndex];
1195 cross *= (currentModel->CrossSectionPerVo 1627 cross *= (currentModel->CrossSectionPerVolume(currentMaterial, particle, kineticEnergy,
1196 1628 (*theCuts)[currentCoupleIndex]));
1197 } 1629 }
1198 return std::max(cross, 0.0); 1630 return std::max(cross, 0.0);
1199 } 1631 }
1200 1632
1201 //....oooOO0OOooo........oooOO0OOooo........o 1633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1202 1634
1203 G4double G4VEnergyLossProcess::MeanFreePath(c 1635 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
1204 { 1636 {
1205 DefineMaterial(track.GetMaterialCutsCouple( 1637 DefineMaterial(track.GetMaterialCutsCouple());
1206 const G4double kinEnergy = track.GetKine 1638 const G4double kinEnergy = track.GetKineticEnergy();
1207 const G4double logKinEnergy = track.GetDyna 1639 const G4double logKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
1208 const G4double cs = GetLambdaForScaledEnerg 1640 const G4double cs = GetLambdaForScaledEnergy(kinEnergy * massRatio,
1209 1641 logKinEnergy + logMassRatio);
1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1642 return (0.0 < cs) ? 1.0/cs : DBL_MAX;
1211 } 1643 }
1212 1644
1213 //....oooOO0OOooo........oooOO0OOooo........o 1645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1214 1646
1215 G4double G4VEnergyLossProcess::ContinuousStep 1647 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
1216 1648 G4double x, G4double y,
1217 1649 G4double& z)
1218 { 1650 {
1219 return AlongStepGetPhysicalInteractionLengt 1651 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &aGPILSelection);
1220 } 1652 }
1221 1653
1222 //....oooOO0OOooo........oooOO0OOooo........o 1654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1223 1655
1224 G4double G4VEnergyLossProcess::GetMeanFreePat 1656 G4double G4VEnergyLossProcess::GetMeanFreePath(
1225 const G4Track& t 1657 const G4Track& track,
1226 G4double, 1658 G4double,
1227 G4ForceCondition 1659 G4ForceCondition* condition)
1228 1660
1229 { 1661 {
1230 *condition = NotForced; 1662 *condition = NotForced;
1231 return MeanFreePath(track); 1663 return MeanFreePath(track);
1232 } 1664 }
1233 1665
1234 //....oooOO0OOooo........oooOO0OOooo........o 1666 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1235 1667
1236 G4double G4VEnergyLossProcess::GetContinuousS 1668 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
1237 const G4Track&, 1669 const G4Track&,
1238 G4double, G4double, G4double& 1670 G4double, G4double, G4double&)
1239 { 1671 {
1240 return DBL_MAX; 1672 return DBL_MAX;
1241 } 1673 }
1242 1674
1243 //....oooOO0OOooo........oooOO0OOooo........o 1675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1244 1676
1245 G4PhysicsVector* 1677 G4PhysicsVector*
1246 G4VEnergyLossProcess::LambdaPhysicsVector(con 1678 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple,
1247 G4d 1679 G4double)
1248 { 1680 {
1249 DefineMaterial(couple); 1681 DefineMaterial(couple);
1250 G4PhysicsVector* v = (*theLambdaTable)[base 1682 G4PhysicsVector* v = (*theLambdaTable)[basedCoupleIndex];
1251 return new G4PhysicsVector(*v); 1683 return new G4PhysicsVector(*v);
1252 } 1684 }
1253 1685
1254 //....oooOO0OOooo........oooOO0OOooo........o 1686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1255 1687
1256 void 1688 void
1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT 1689 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
1258 { 1690 {
1259 if(1 < verboseLevel) { <<
1260 G4cout << "### Set DEDX table " << p << " <<
1261 << " " << theDEDXunRestrictedTable << <<
1262 << " for " << particle->GetParticl <<
1263 << " and process " << GetProcessNa <<
1264 << " type=" << tType << " isIonisation:" <<
1265 } <<
1266 if(fTotal == tType) { 1691 if(fTotal == tType) {
1267 theDEDXunRestrictedTable = p; 1692 theDEDXunRestrictedTable = p;
>> 1693
1268 } else if(fRestricted == tType) { 1694 } else if(fRestricted == tType) {
>> 1695 /*
>> 1696 G4cout<< "G4VEnergyLossProcess::SetDEDXTable "
>> 1697 << particle->GetParticleName()
>> 1698 << " oldTable " << theDEDXTable << " newTable " << p
>> 1699 << " ion " << theIonisationTable
>> 1700 << " IsMaster " << isMaster
>> 1701 << " " << GetProcessName() << G4endl;
>> 1702 G4cout << (*p) << G4endl;
>> 1703 */
1269 theDEDXTable = p; 1704 theDEDXTable = p;
1270 if(isMaster && nullptr == baseParticle) { <<
1271 theData->UpdateTable(theDEDXTable, 0); <<
1272 } <<
1273 } else if(fIsIonisation == tType) { 1705 } else if(fIsIonisation == tType) {
>> 1706 /*
>> 1707 G4cout<< "G4VEnergyLossProcess::SetIonisationTable "
>> 1708 << particle->GetParticleName()
>> 1709 << " oldTable " << theDEDXTable << " newTable " << p
>> 1710 << " ion " << theIonisationTable
>> 1711 << " IsMaster " << isMaster
>> 1712 << " " << GetProcessName() << G4endl;
>> 1713 */
1274 theIonisationTable = p; 1714 theIonisationTable = p;
1275 if(isMaster && nullptr == baseParticle) { <<
1276 theData->UpdateTable(theIonisationTable <<
1277 } <<
1278 } 1715 }
1279 } 1716 }
1280 1717
1281 //....oooOO0OOooo........oooOO0OOooo........o 1718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1282 1719
1283 void G4VEnergyLossProcess::SetCSDARangeTable( 1720 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
1284 { 1721 {
1285 theCSDARangeTable = p; << 1722 theCSDARangeTable = p;
>> 1723 if(1 < verboseLevel) {
>> 1724 G4cout << "### Set CSDA Range table " << p
>> 1725 << " for " << particle->GetParticleName()
>> 1726 << " and process " << GetProcessName() << G4endl;
>> 1727 }
1286 } 1728 }
1287 1729
1288 //....oooOO0OOooo........oooOO0OOooo........o 1730 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1289 1731
1290 void G4VEnergyLossProcess::SetRangeTableForLo 1732 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
1291 { 1733 {
1292 theRangeTableForLoss = p; 1734 theRangeTableForLoss = p;
>> 1735 if(1 < verboseLevel) {
>> 1736 G4cout << "### Set Range table " << p
>> 1737 << " for " << particle->GetParticleName()
>> 1738 << " and process " << GetProcessName() << G4endl;
>> 1739 }
>> 1740 }
>> 1741
>> 1742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 1743
>> 1744 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p)
>> 1745 {
>> 1746 theSecondaryRangeTable = p;
>> 1747 if(1 < verboseLevel) {
>> 1748 G4cout << "### Set SecondaryRange table " << p
>> 1749 << " for " << particle->GetParticleName()
>> 1750 << " and process " << GetProcessName() << G4endl;
>> 1751 }
1293 } 1752 }
1294 1753
1295 //....oooOO0OOooo........oooOO0OOooo........o 1754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1296 1755
1297 void G4VEnergyLossProcess::SetInverseRangeTab 1756 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
1298 { 1757 {
1299 theInverseRangeTable = p; 1758 theInverseRangeTable = p;
>> 1759 if(1 < verboseLevel) {
>> 1760 G4cout << "### Set InverseRange table " << p
>> 1761 << " for " << particle->GetParticleName()
>> 1762 << " and process " << GetProcessName() << G4endl;
>> 1763 }
1300 } 1764 }
1301 1765
1302 //....oooOO0OOooo........oooOO0OOooo........o 1766 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1303 1767
1304 void G4VEnergyLossProcess::SetLambdaTable(G4P 1768 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
1305 { 1769 {
1306 if(1 < verboseLevel) { 1770 if(1 < verboseLevel) {
1307 G4cout << "### Set Lambda table " << p << << 1771 G4cout << "### Set Lambda table " << p
1308 << " for " << particle->GetParticl 1772 << " for " << particle->GetParticleName()
1309 << " and process " << GetProcessNa 1773 << " and process " << GetProcessName() << G4endl;
>> 1774 //G4cout << *p << G4endl;
1310 } 1775 }
1311 theLambdaTable = p; << 1776 theLambdaTable = p;
1312 tablesAreBuilt = true; 1777 tablesAreBuilt = true;
1313 1778
1314 if(isMaster && nullptr != p) { << 1779 G4LossTableBuilder* bld = lManager->GetTableBuilder();
1315 delete theEnergyOfCrossSectionMax; << 1780 theDensityFactor = bld->GetDensityFactors();
1316 theEnergyOfCrossSectionMax = nullptr; << 1781 theDensityIdx = bld->GetCoupleIndexes();
1317 if(fEmTwoPeaks == fXSType) { << 1782
1318 if(nullptr != fXSpeaks) { << 1783 if(isMaster && nullptr == baseParticle &&
1319 for(auto & ptr : *fXSpeaks) { delete ptr; } << 1784 nullptr != theLambdaTable && fEmTwoPeaks == fXSType) {
1320 delete fXSpeaks; << 1785
>> 1786 size_t n = theLambdaTable->length();
>> 1787
>> 1788 G4double e, ss, xs, ee, e1peak, xs1peak, e1deep, e2peak, e2deep, xs2peak;
>> 1789
>> 1790 // first loop on existing vectors
>> 1791 for (size_t i=0; i<n; ++i) {
>> 1792 const G4PhysicsVector* pv = (*theLambdaTable)[i];
>> 1793 ee = xs = xs1peak = xs2peak = 0.0;
>> 1794 e1peak = e1deep = e2peak = e2deep = DBL_MAX;
>> 1795 if(nullptr != pv) {
>> 1796 size_t nb = pv->GetVectorLength();
>> 1797 for (size_t j=0; j<nb; ++j) {
>> 1798 e = pv->Energy(j);
>> 1799 ss = (*pv)(j);
>> 1800 // find out 1st peak
>> 1801 if(e1peak == DBL_MAX) {
>> 1802 if(ss >= xs) {
>> 1803 xs = ss;
>> 1804 ee = e;
>> 1805 continue;
>> 1806 } else {
>> 1807 e1peak = ee;
>> 1808 xs1peak = xs;
>> 1809 }
>> 1810 }
>> 1811 // find out the deep
>> 1812 if(e1deep == DBL_MAX) {
>> 1813 if(ss <= xs) {
>> 1814 xs = ss;
>> 1815 ee = e;
>> 1816 continue;
>> 1817 } else {
>> 1818 e1deep = ee;
>> 1819 }
>> 1820 }
>> 1821 // find out 2nd peak
>> 1822 if(e2peak == DBL_MAX) {
>> 1823 if(ss >= xs) {
>> 1824 xs = ss;
>> 1825 ee = e;
>> 1826 continue;
>> 1827 } else {
>> 1828 e2peak = ee;
>> 1829 xs2peak = xs;
>> 1830 }
>> 1831 }
>> 1832 if(e2deep == DBL_MAX) {
>> 1833 if(ss <= xs) {
>> 1834 xs = ss;
>> 1835 ee = e;
>> 1836 continue;
>> 1837 } else {
>> 1838 e2deep = ee;
>> 1839 break;
>> 1840 }
>> 1841 }
>> 1842 }
>> 1843 }
>> 1844 G4TwoPeaksXS* x = (*fXSpeaks)[i];
>> 1845 if(nullptr == x) {
>> 1846 x = new G4TwoPeaksXS();
>> 1847 (*fXSpeaks)[i] = x;
>> 1848 }
>> 1849 x->e1peak = e1peak;
>> 1850 x->e1deep = e1deep;
>> 1851 x->e2peak = e2peak;
>> 1852 x->e2deep = e2deep;
>> 1853
>> 1854 if(1 < verboseLevel) {
>> 1855 G4cout << "For " << particle->GetParticleName()
>> 1856 << " index= " << i << " data:\n" << " E1peak=" << e1peak
>> 1857 << " xs1= " << xs1peak << " E1deep=" << e1deep
>> 1858 << " E2peak=" << e2peak << " xs2=" << xs2peak
>> 1859 << " E2deep=" << e2deep << G4endl;
1321 } 1860 }
1322 G4LossTableBuilder* bld = lManager->Get <<
1323 fXSpeaks = G4EmUtility::FillPeaksStruct <<
1324 if(nullptr == fXSpeaks) { fXSType = fEm <<
1325 } 1861 }
1326 if(fXSType == fEmOnePeak) { << 1862 // second loop using base materials
1327 theEnergyOfCrossSectionMax = G4EmUtilit << 1863 for (size_t i=0; i<n; ++i) {
1328 if(nullptr == theEnergyOfCrossSectionMa << 1864 const G4PhysicsVector* pv = (*theLambdaTable)[i];
>> 1865 if (nullptr == pv) {
>> 1866 G4int j = (*theDensityIdx)[i];
>> 1867 G4TwoPeaksXS* x = (*fXSpeaks)[i];
>> 1868 G4TwoPeaksXS* y = (*fXSpeaks)[j];
>> 1869 if(nullptr == x) {
>> 1870 x = new G4TwoPeaksXS();
>> 1871 (*fXSpeaks)[i] = x;
>> 1872 }
>> 1873 x->e1peak = y->e1peak;
>> 1874 x->e1deep = y->e1deep;
>> 1875 x->e2peak = y->e2peak;
>> 1876 x->e2deep = y->e2deep;
>> 1877 }
1329 } 1878 }
1330 } 1879 }
1331 } 1880 }
1332 1881
1333 //....oooOO0OOooo........oooOO0OOooo........o 1882 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1334 1883
1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe <<
1336 { <<
1337 theEnergyOfCrossSectionMax = p; <<
1338 } <<
1339 <<
1340 //....oooOO0OOooo........oooOO0OOooo........o <<
1341 <<
1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: 1884 void G4VEnergyLossProcess::SetTwoPeaksXS(std::vector<G4TwoPeaksXS*>* ptr)
1343 { 1885 {
1344 fXSpeaks = ptr; 1886 fXSpeaks = ptr;
1345 } 1887 }
1346 1888
1347 //....oooOO0OOooo........oooOO0OOooo........o 1889 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1348 1890
1349 const G4Element* G4VEnergyLossProcess::GetCur 1891 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
1350 { 1892 {
1351 return (nullptr != currentModel) << 1893 return (nullptr != currentModel) ? currentModel->GetCurrentElement() : nullptr;
1352 ? currentModel->GetCurrentElement(current <<
1353 } 1894 }
1354 1895
1355 //....oooOO0OOooo........oooOO0OOooo........o 1896 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1356 1897
1357 void G4VEnergyLossProcess::SetCrossSectionBia 1898 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f,
1358 1899 G4bool flag)
1359 { 1900 {
1360 if(f > 0.0) { 1901 if(f > 0.0) {
1361 biasFactor = f; 1902 biasFactor = f;
1362 weightFlag = flag; 1903 weightFlag = flag;
1363 if(1 < verboseLevel) { 1904 if(1 < verboseLevel) {
1364 G4cout << "### SetCrossSectionBiasingFa 1905 G4cout << "### SetCrossSectionBiasingFactor: for "
1365 << " process " << GetProcessName 1906 << " process " << GetProcessName()
1366 << " biasFactor= " << f << " wei 1907 << " biasFactor= " << f << " weightFlag= " << flag
1367 << G4endl; 1908 << G4endl;
1368 } 1909 }
1369 } 1910 }
1370 } 1911 }
1371 1912
1372 //....oooOO0OOooo........oooOO0OOooo........o 1913 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1373 1914
1374 void G4VEnergyLossProcess::ActivateForcedInte 1915 void G4VEnergyLossProcess::ActivateForcedInteraction(G4double length,
1375 1916 const G4String& region,
1376 1917 G4bool flag)
1377 { 1918 {
1378 if(nullptr == biasManager) { biasManager = 1919 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1379 if(1 < verboseLevel) { 1920 if(1 < verboseLevel) {
1380 G4cout << "### ActivateForcedInteraction: 1921 G4cout << "### ActivateForcedInteraction: for "
1381 << " process " << GetProcessName() 1922 << " process " << GetProcessName()
1382 << " length(mm)= " << length/mm 1923 << " length(mm)= " << length/mm
1383 << " in G4Region <" << region 1924 << " in G4Region <" << region
1384 << "> weightFlag= " << flag 1925 << "> weightFlag= " << flag
1385 << G4endl; 1926 << G4endl;
1386 } 1927 }
1387 weightFlag = flag; 1928 weightFlag = flag;
1388 biasManager->ActivateForcedInteraction(leng 1929 biasManager->ActivateForcedInteraction(length, region);
1389 } 1930 }
1390 1931
1391 //....oooOO0OOooo........oooOO0OOooo........o 1932 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1392 1933
1393 void 1934 void
1394 G4VEnergyLossProcess::ActivateSecondaryBiasin 1935 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region,
1395 1936 G4double factor,
1396 1937 G4double energyLimit)
1397 { 1938 {
1398 if (0.0 <= factor) { 1939 if (0.0 <= factor) {
1399 // Range cut can be applied only for e- 1940 // Range cut can be applied only for e-
1400 if(0.0 == factor && secondaryParticle != 1941 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
1401 { return; } 1942 { return; }
1402 1943
1403 if(nullptr == biasManager) { biasManager 1944 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1404 biasManager->ActivateSecondaryBiasing(reg 1945 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
1405 if(1 < verboseLevel) { 1946 if(1 < verboseLevel) {
1406 G4cout << "### ActivateSecondaryBiasing 1947 G4cout << "### ActivateSecondaryBiasing: for "
1407 << " process " << GetProcessName 1948 << " process " << GetProcessName()
1408 << " factor= " << factor 1949 << " factor= " << factor
1409 << " in G4Region <" << region 1950 << " in G4Region <" << region
1410 << "> energyLimit(MeV)= " << ene 1951 << "> energyLimit(MeV)= " << energyLimit/MeV
1411 << G4endl; 1952 << G4endl;
1412 } 1953 }
1413 } 1954 }
1414 } 1955 }
1415 1956
1416 //....oooOO0OOooo........oooOO0OOooo........o 1957 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1417 1958
1418 void G4VEnergyLossProcess::SetIonisation(G4bo 1959 void G4VEnergyLossProcess::SetIonisation(G4bool val)
1419 { 1960 {
1420 isIonisation = val; 1961 isIonisation = val;
1421 aGPILSelection = (val) ? CandidateForSelect 1962 aGPILSelection = (val) ? CandidateForSelection : NotCandidateForSelection;
1422 } 1963 }
1423 1964
1424 //....oooOO0OOooo........oooOO0OOooo........o 1965 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1425 1966
1426 void G4VEnergyLossProcess::SetLinearLossLimi 1967 void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
1427 { 1968 {
1428 if(0.0 < val && val < 1.0) { 1969 if(0.0 < val && val < 1.0) {
1429 linLossLimit = val; 1970 linLossLimit = val;
1430 actLinLossLimit = true; 1971 actLinLossLimit = true;
1431 } else { PrintWarning("SetLinearLossLimit", 1972 } else { PrintWarning("SetLinearLossLimit", val); }
1432 } 1973 }
1433 1974
1434 //....oooOO0OOooo........oooOO0OOooo........o 1975 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1435 1976
1436 void G4VEnergyLossProcess::SetStepFunction(G4 1977 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
1437 { 1978 {
1438 if(0.0 < v1 && 0.0 < v2) { 1979 if(0.0 < v1 && 0.0 < v2) {
1439 dRoverRange = std::min(1.0, v1); 1980 dRoverRange = std::min(1.0, v1);
1440 finalRange = std::min(v2, 1.e+50); 1981 finalRange = std::min(v2, 1.e+50);
1441 } else { 1982 } else {
1442 PrintWarning("SetStepFunctionV1", v1); 1983 PrintWarning("SetStepFunctionV1", v1);
1443 PrintWarning("SetStepFunctionV2", v2); 1984 PrintWarning("SetStepFunctionV2", v2);
1444 } 1985 }
1445 } 1986 }
1446 1987
1447 //....oooOO0OOooo........oooOO0OOooo........o 1988 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1448 1989
1449 void G4VEnergyLossProcess::SetLowestEnergyLim 1990 void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val)
1450 { 1991 {
1451 if(1.e-18 < val && val < 1.e+50) { lowestKi 1992 if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; }
1452 else { PrintWarning("SetLowestEnergyLimit", 1993 else { PrintWarning("SetLowestEnergyLimit", val); }
1453 } 1994 }
1454 1995
1455 //....oooOO0OOooo........oooOO0OOooo........o 1996 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1456 1997
1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i 1998 void G4VEnergyLossProcess::SetDEDXBinning(G4int n)
1458 { 1999 {
1459 if(2 < n && n < 1000000000) { 2000 if(2 < n && n < 1000000000) {
1460 nBins = n; 2001 nBins = n;
1461 actBinning = true; 2002 actBinning = true;
1462 } else { 2003 } else {
1463 G4double e = (G4double)n; 2004 G4double e = (G4double)n;
1464 PrintWarning("SetDEDXBinning", e); 2005 PrintWarning("SetDEDXBinning", e);
1465 } 2006 }
1466 } 2007 }
1467 2008
1468 //....oooOO0OOooo........oooOO0OOooo........o 2009 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1469 2010
1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 2011 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
1471 { 2012 {
1472 if(1.e-18 < e && e < maxKinEnergy) { 2013 if(1.e-18 < e && e < maxKinEnergy) {
1473 minKinEnergy = e; 2014 minKinEnergy = e;
1474 actMinKinEnergy = true; 2015 actMinKinEnergy = true;
1475 } else { PrintWarning("SetMinKinEnergy", e) 2016 } else { PrintWarning("SetMinKinEnergy", e); }
1476 } 2017 }
1477 2018
1478 //....oooOO0OOooo........oooOO0OOooo........o 2019 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1479 2020
1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 2021 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
1481 { 2022 {
1482 if(minKinEnergy < e && e < 1.e+50) { 2023 if(minKinEnergy < e && e < 1.e+50) {
1483 maxKinEnergy = e; 2024 maxKinEnergy = e;
1484 actMaxKinEnergy = true; 2025 actMaxKinEnergy = true;
1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS 2026 if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; }
1486 } else { PrintWarning("SetMaxKinEnergy", e) 2027 } else { PrintWarning("SetMaxKinEnergy", e); }
1487 } 2028 }
1488 2029
1489 //....oooOO0OOooo........oooOO0OOooo........o 2030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1490 2031
1491 void G4VEnergyLossProcess::PrintWarning(const 2032 void G4VEnergyLossProcess::PrintWarning(const G4String& tit, G4double val) const
1492 { 2033 {
1493 G4String ss = "G4VEnergyLossProcess::" + ti 2034 G4String ss = "G4VEnergyLossProcess::" + tit;
1494 G4ExceptionDescription ed; 2035 G4ExceptionDescription ed;
1495 ed << "Parameter is out of range: " << val 2036 ed << "Parameter is out of range: " << val
1496 << " it will have no effect!\n" << " Pr 2037 << " it will have no effect!\n" << " Process "
1497 << GetProcessName() << " nbins= " << nB 2038 << GetProcessName() << " nbins= " << nBins
1498 << " Emin(keV)= " << minKinEnergy/keV 2039 << " Emin(keV)= " << minKinEnergy/keV
1499 << " Emax(GeV)= " << maxKinEnergy/GeV; 2040 << " Emax(GeV)= " << maxKinEnergy/GeV;
1500 G4Exception(ss, "em0044", JustWarning, ed); 2041 G4Exception(ss, "em0044", JustWarning, ed);
1501 } 2042 }
1502 2043
1503 //....oooOO0OOooo........oooOO0OOooo........o 2044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1504 2045
1505 void G4VEnergyLossProcess::ProcessDescription 2046 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const
1506 { 2047 {
1507 if(nullptr != particle) { StreamInfo(out, * 2048 if(nullptr != particle) { StreamInfo(out, *particle, true); }
1508 } 2049 }
1509 2050
1510 //....oooOO0OOooo........oooOO0OOooo........o 2051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1511 2052