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