| 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: G4VEmProcess 31 // File name: G4VEmProcess
32 // 32 //
33 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 // 34 //
35 // Creation date: 01.10.2003 35 // Creation date: 01.10.2003
36 // 36 //
37 // Modifications: by V.Ivanchenko 37 // Modifications: by V.Ivanchenko
38 // 38 //
39 // Class Description: based class for discrete 39 // Class Description: based class for discrete and rest/discrete EM processes
40 // 40 //
41 41
42 // ------------------------------------------- 42 // -------------------------------------------------------------------
43 // 43 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 46
47 #include "G4VEmProcess.hh" 47 #include "G4VEmProcess.hh"
48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
50 #include "G4ProcessManager.hh" 50 #include "G4ProcessManager.hh"
51 #include "G4LossTableManager.hh" 51 #include "G4LossTableManager.hh"
52 #include "G4LossTableBuilder.hh" 52 #include "G4LossTableBuilder.hh"
53 #include "G4Step.hh" 53 #include "G4Step.hh"
54 #include "G4ParticleDefinition.hh" 54 #include "G4ParticleDefinition.hh"
55 #include "G4VEmModel.hh" 55 #include "G4VEmModel.hh"
56 #include "G4DataVector.hh" 56 #include "G4DataVector.hh"
57 #include "G4PhysicsTable.hh" 57 #include "G4PhysicsTable.hh"
58 #include "G4EmDataHandler.hh" 58 #include "G4EmDataHandler.hh"
59 #include "G4PhysicsLogVector.hh" 59 #include "G4PhysicsLogVector.hh"
60 #include "G4VParticleChange.hh" 60 #include "G4VParticleChange.hh"
61 #include "G4ProductionCutsTable.hh" 61 #include "G4ProductionCutsTable.hh"
62 #include "G4Region.hh" 62 #include "G4Region.hh"
63 #include "G4Gamma.hh" 63 #include "G4Gamma.hh"
64 #include "G4Electron.hh" 64 #include "G4Electron.hh"
65 #include "G4Positron.hh" 65 #include "G4Positron.hh"
66 #include "G4PhysicsTableHelper.hh" 66 #include "G4PhysicsTableHelper.hh"
67 #include "G4EmBiasingManager.hh" 67 #include "G4EmBiasingManager.hh"
68 #include "G4EmParameters.hh" <<
69 #include "G4EmProcessSubType.hh" <<
70 #include "G4EmTableUtil.hh" <<
71 #include "G4EmUtility.hh" <<
72 #include "G4DNAModelSubType.hh" <<
73 #include "G4GenericIon.hh" 68 #include "G4GenericIon.hh"
74 #include "G4Log.hh" 69 #include "G4Log.hh"
75 #include <iostream> 70 #include <iostream>
76 71
77 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 73
79 G4VEmProcess::G4VEmProcess(const G4String& nam 74 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
80 G4VDiscreteProcess(name, type) << 75 G4VDiscreteProcess(name, type),
>> 76 secondaryParticle(nullptr),
>> 77 buildLambdaTable(true),
>> 78 numberOfModels(0),
>> 79 theLambdaTable(nullptr),
>> 80 theLambdaTablePrim(nullptr),
>> 81 integral(false),
>> 82 applyCuts(false),
>> 83 startFromNull(false),
>> 84 splineFlag(true),
>> 85 isIon(false),
>> 86 currentCouple(nullptr),
>> 87 isTheMaster(true),
>> 88 masterProc(nullptr),
>> 89 theData(nullptr),
>> 90 currentModel(nullptr),
>> 91 particle(nullptr),
>> 92 currentParticle(nullptr)
81 { 93 {
82 theParameters = G4EmParameters::Instance(); 94 theParameters = G4EmParameters::Instance();
83 SetVerboseLevel(1); 95 SetVerboseLevel(1);
84 96
85 // Size of tables << 97 // Size of tables assuming spline
86 minKinEnergy = 0.1*CLHEP::keV; << 98 minKinEnergy = 0.1*keV;
87 maxKinEnergy = 100.0*CLHEP::TeV; << 99 maxKinEnergy = 100.0*TeV;
>> 100 nLambdaBins = 84;
>> 101 minKinEnergyPrim = DBL_MAX;
>> 102 actBinning = actSpline = actMinKinEnergy = actMaxKinEnergy = false;
88 103
89 // default lambda factor 104 // default lambda factor
90 invLambdaFactor = 1.0/lambdaFactor; << 105 lambdaFactor = 0.8;
>> 106 logLambdaFactor = G4Log(lambdaFactor);
>> 107
>> 108 // default limit on polar angle
>> 109 biasFactor = fFactor = 1.0;
91 110
92 // particle types 111 // particle types
93 theGamma = G4Gamma::Gamma(); << 112 theGamma = G4Gamma::Gamma();
94 theElectron = G4Electron::Electron(); << 113 theElectron = G4Electron::Electron();
95 thePositron = G4Positron::Positron(); << 114 thePositron = G4Positron::Positron();
>> 115
>> 116 theCuts = theCutsGamma = theCutsElectron = theCutsPositron = nullptr;
96 117
97 pParticleChange = &fParticleChange; 118 pParticleChange = &fParticleChange;
98 fParticleChange.SetSecondaryWeightByProcess( 119 fParticleChange.SetSecondaryWeightByProcess(true);
99 secParticles.reserve(5); 120 secParticles.reserve(5);
100 121
>> 122 baseMaterial = currentMaterial = nullptr;
>> 123
>> 124 preStepLambda = preStepKinEnergy = 0.0;
>> 125 preStepLogKinEnergy = LOG_EKIN_MIN;
>> 126 mfpKinEnergy = DBL_MAX;
>> 127 massRatio = 1.0;
>> 128
>> 129 currentCoupleIndex = basedCoupleIndex = 0;
>> 130
101 modelManager = new G4EmModelManager(); 131 modelManager = new G4EmModelManager();
>> 132 biasManager = nullptr;
>> 133 biasFlag = false;
>> 134 weightFlag = false;
102 lManager = G4LossTableManager::Instance(); 135 lManager = G4LossTableManager::Instance();
103 lManager->Register(this); 136 lManager->Register(this);
104 isTheMaster = lManager->IsMaster(); <<
105 G4LossTableBuilder* bld = lManager->GetTable 137 G4LossTableBuilder* bld = lManager->GetTableBuilder();
106 theDensityFactor = bld->GetDensityFactors(); 138 theDensityFactor = bld->GetDensityFactors();
107 theDensityIdx = bld->GetCoupleIndexes(); 139 theDensityIdx = bld->GetCoupleIndexes();
>> 140
>> 141 secID = fluoID = augerID = biasID = -1;
>> 142 mainSecondaries = 100;
>> 143 if("phot" == GetProcessName() || "compt" == GetProcessName()
>> 144 || "e-_G4DNAIonisation" == GetProcessName()
>> 145 || "hydrogen_G4DNAIonisation" == GetProcessName()
>> 146 || "helium_G4DNAIonisation" == GetProcessName()
>> 147 || "alpha_G4DNAIonisation" == GetProcessName()
>> 148 || "alpha+_G4DNAIonisation" == GetProcessName()
>> 149 || "proton_G4DNAIonisation" == GetProcessName()
>> 150 || "GenericIon_G4DNAIonisation" == GetProcessName() )
>> 151 {
>> 152 mainSecondaries = 1;
>> 153 }
108 } 154 }
109 155
110 //....oooOO0OOooo........oooOO0OOooo........oo 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 157
112 G4VEmProcess::~G4VEmProcess() 158 G4VEmProcess::~G4VEmProcess()
113 { 159 {
>> 160 /*
>> 161 if(1 < verboseLevel) {
>> 162 G4cout << "G4VEmProcess destruct " << GetProcessName()
>> 163 << " " << this << " " << theLambdaTable <<G4endl;
>> 164 }
>> 165 */
114 if(isTheMaster) { 166 if(isTheMaster) {
115 delete theData; 167 delete theData;
116 delete theEnergyOfCrossSectionMax; << 168 theData = nullptr;
117 } 169 }
118 delete modelManager; 170 delete modelManager;
119 delete biasManager; 171 delete biasManager;
120 lManager->DeRegister(this); 172 lManager->DeRegister(this);
121 } 173 }
122 174
123 //....oooOO0OOooo........oooOO0OOooo........oo 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 176
125 void G4VEmProcess::AddEmModel(G4int order, G4V << 177 void G4VEmProcess::Clear()
>> 178 {
>> 179 currentCouple = nullptr;
>> 180 preStepLambda = 0.0;
>> 181 }
>> 182
>> 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 184
>> 185 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
>> 186 const G4Material*)
>> 187 {
>> 188 return 0.0;
>> 189 }
>> 190
>> 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 192
>> 193 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
126 const G4Region* 194 const G4Region* region)
127 { 195 {
128 if(nullptr == ptr) { return; } <<
129 G4VEmFluctuationModel* fm = nullptr; 196 G4VEmFluctuationModel* fm = nullptr;
130 modelManager->AddEmModel(order, ptr, fm, reg << 197 modelManager->AddEmModel(order, p, fm, region);
131 ptr->SetParticleChange(pParticleChange); << 198 if(p) { p->SetParticleChange(pParticleChange); }
132 } 199 }
133 200
134 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 202
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, 203 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, G4int)
137 { 204 {
138 if(nullptr == ptr) { return; } << 205 for(auto & em : emModels) { if(em == ptr) { return; } }
139 if(!emModels.empty()) { << 206 emModels.push_back(ptr);
140 for(auto & em : emModels) { if(em == ptr) << 207 }
141 } << 208
142 emModels.push_back(ptr); << 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 210
>> 211 G4VEmModel* G4VEmProcess::EmModel(size_t index) const
>> 212 {
>> 213 return (index < emModels.size()) ? emModels[index] : nullptr;
>> 214 }
>> 215
>> 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 217
>> 218 void G4VEmProcess::UpdateEmModel(const G4String& nam,
>> 219 G4double emin, G4double emax)
>> 220 {
>> 221 modelManager->UpdateEmModel(nam, emin, emax);
>> 222 }
>> 223
>> 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 225
>> 226 G4int G4VEmProcess::GetNumberOfModels() const
>> 227 {
>> 228 return modelManager->NumberOfModels();
>> 229 }
>> 230
>> 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 232
>> 233 G4int G4VEmProcess::GetNumberOfRegionModels(size_t couple_index) const
>> 234 {
>> 235 return modelManager->NumberOfRegionModels(couple_index);
>> 236 }
>> 237
>> 238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 239
>> 240 G4VEmModel* G4VEmProcess::GetRegionModel(G4int idx, size_t couple_index) const
>> 241 {
>> 242 return modelManager->GetRegionModel(idx, couple_index);
>> 243 }
>> 244
>> 245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 246
>> 247 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver) const
>> 248 {
>> 249 return modelManager->GetModel(idx, ver);
143 } 250 }
144 251
145 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 253
147 void G4VEmProcess::PreparePhysicsTable(const G 254 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
148 { 255 {
149 if(nullptr == particle) { SetParticle(&part) << 256 isTheMaster = lManager->IsMaster();
>> 257
>> 258 if(!particle) { SetParticle(&part); }
150 259
151 if(part.GetParticleType() == "nucleus" && 260 if(part.GetParticleType() == "nucleus" &&
152 part.GetParticleSubType() == "generic") { 261 part.GetParticleSubType() == "generic") {
153 262
154 G4String pname = part.GetParticleName(); 263 G4String pname = part.GetParticleName();
155 if(pname != "deuteron" && pname != "triton 264 if(pname != "deuteron" && pname != "triton" &&
156 pname != "He3" && pname != "alpha" && p << 265 pname != "alpha" && pname != "He3" &&
157 pname != "helium" && pname != "hydrogen << 266 pname != "alpha+" && pname != "helium" &&
>> 267 pname != "hydrogen") {
158 268
159 particle = G4GenericIon::GenericIon(); 269 particle = G4GenericIon::GenericIon();
160 isIon = true; 270 isIon = true;
161 } 271 }
162 } 272 }
>> 273
>> 274 if(1 < verboseLevel) {
>> 275 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
>> 276 << GetProcessName()
>> 277 << " and particle " << part.GetParticleName()
>> 278 << " local particle " << particle->GetParticleName()
>> 279 << G4endl;
>> 280 }
>> 281
163 if(particle != &part) { return; } 282 if(particle != &part) { return; }
164 283
165 lManager->PreparePhysicsTable(&part, this); << 284 G4LossTableBuilder* bld = lManager->GetTableBuilder();
166 285
167 // for new run << 286 lManager->PreparePhysicsTable(&part, this, isTheMaster);
168 currentCouple = nullptr; <<
169 preStepLambda = 0.0; <<
170 fLambdaEnergy = 0.0; <<
171 287
>> 288 Clear();
172 InitialiseProcess(particle); 289 InitialiseProcess(particle);
173 290
174 G4LossTableBuilder* bld = lManager->GetTable <<
175 const G4ProductionCutsTable* theCoupleTable= 291 const G4ProductionCutsTable* theCoupleTable=
176 G4ProductionCutsTable::GetProductionCutsTa 292 G4ProductionCutsTable::GetProductionCutsTable();
177 theCutsGamma = theCoupleTable->GetEnergyC << 293 size_t n = theCoupleTable->GetTableSize();
178 theCutsElectron = theCoupleTable->GetEnergyC << 294
179 theCutsPositron = theCoupleTable->GetEnergyC << 295 theEnergyOfCrossSectionMax.resize(n, 0.0);
>> 296 theCrossSectionMax.resize(n, DBL_MAX);
180 297
181 // initialisation of the process 298 // initialisation of the process
182 if(!actMinKinEnergy) { minKinEnergy = thePar 299 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
183 if(!actMaxKinEnergy) { maxKinEnergy = thePar 300 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
>> 301 if(!actSpline) { splineFlag = theParameters->Spline(); }
184 302
>> 303 if(isTheMaster) {
>> 304 SetVerboseLevel(theParameters->Verbose());
>> 305 if(!theData) { theData = new G4EmDataHandler(2); }
>> 306 } else {
>> 307 SetVerboseLevel(theParameters->WorkerVerbose());
>> 308 }
185 applyCuts = theParameters->ApplyCuts() 309 applyCuts = theParameters->ApplyCuts();
186 lambdaFactor = theParameters->LambdaFacto 310 lambdaFactor = theParameters->LambdaFactor();
187 invLambdaFactor = 1.0/lambdaFactor; << 311 logLambdaFactor = G4Log(lambdaFactor);
188 theParameters->DefineRegParamForEM(this); 312 theParameters->DefineRegParamForEM(this);
189 313
190 // integral option may be disabled << 314 // initialisation of models
191 if(!theParameters->Integral()) { fXSType = f << 315 numberOfModels = modelManager->NumberOfModels();
192 << 316 for(G4int i=0; i<numberOfModels; ++i) {
193 // prepare tables << 317 G4VEmModel* mod = modelManager->GetModel(i);
194 if(isTheMaster) { << 318 if(0 == i) { currentModel = mod; }
195 if(nullptr == theData) { theData = new G4E << 319 mod->SetPolarAngleLimit(theParameters->MscThetaLimit());
196 << 320 mod->SetMasterThread(isTheMaster);
197 if(buildLambdaTable) { << 321 if(mod->HighEnergyLimit() > maxKinEnergy) {
198 theLambdaTable = theData->MakeTable(0); << 322 mod->SetHighEnergyLimit(maxKinEnergy);
199 bld->InitialiseBaseMaterials(theLambdaTa <<
200 } <<
201 // high energy table <<
202 if(minKinEnergyPrim < maxKinEnergy) { <<
203 theLambdaTablePrim = theData->MakeTable( <<
204 bld->InitialiseBaseMaterials(theLambdaTa <<
205 } 323 }
206 } 324 }
207 // models << 325
208 baseMat = bld->GetBaseMaterialFlag(); << 326 if(lManager->AtomDeexcitation()) { modelManager->SetFluoFlag(true); }
209 numberOfModels = modelManager->NumberOfModel << 327 theCuts = modelManager->Initialise(particle,secondaryParticle,
210 currentModel = modelManager->GetModel(0); << 328 2.,verboseLevel);
211 if(nullptr != lManager->AtomDeexcitation()) << 329 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
212 modelManager->SetFluoFlag(true); << 330 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
>> 331 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
>> 332
>> 333 // prepare tables
>> 334 if(buildLambdaTable && isTheMaster){
>> 335 theLambdaTable = theData->MakeTable(0);
>> 336 bld->InitialiseBaseMaterials(theLambdaTable);
>> 337 }
>> 338 // high energy table
>> 339 if(isTheMaster && minKinEnergyPrim < maxKinEnergy){
>> 340 theLambdaTablePrim = theData->MakeTable(1);
>> 341 bld->InitialiseBaseMaterials(theLambdaTablePrim);
213 } 342 }
>> 343 bld->InitialiseBaseMaterials();
214 // forced biasing 344 // forced biasing
215 if(nullptr != biasManager) { << 345 if(biasManager) {
216 biasManager->Initialise(part, GetProcessNa << 346 biasManager->Initialise(part,GetProcessName(),verboseLevel);
217 biasFlag = false; << 347 biasFlag = false;
218 } << 348 }
219 << 349 // defined ID of secondary particles
220 theCuts = << 350 G4String nam1 = GetProcessName();
221 G4EmTableUtil::PrepareEmProcess(this, part << 351 secID = G4PhysicsModelCatalog::Register(nam1);
222 modelManag << 352 if(100 > mainSecondaries) {
223 secID, tri << 353 G4String nam2 = nam1 + "_fluo" ;
224 verboseLev << 354 G4String nam3 = nam1 + "_auger";
>> 355 G4String nam4 = nam1 + "_split";
>> 356 fluoID = G4PhysicsModelCatalog::Register(nam2);
>> 357 augerID = G4PhysicsModelCatalog::Register(nam3);
>> 358 biasID = G4PhysicsModelCatalog::Register(nam4);
>> 359 }
225 } 360 }
226 361
227 //....oooOO0OOooo........oooOO0OOooo........oo 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 363
229 void G4VEmProcess::BuildPhysicsTable(const G4P 364 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 { 365 {
231 if(nullptr == masterProc) { << 366 if(!masterProc) {
232 if(isTheMaster) { masterProc = this; } 367 if(isTheMaster) { masterProc = this; }
233 else { masterProc = static_cast<const G4VE 368 else { masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess());}
234 } 369 }
235 G4int nModels = modelManager->NumberOfModels << 370
236 G4bool isLocked = theParameters->IsPrintLock << 371 G4String num = part.GetParticleName();
237 G4bool toBuild = (buildLambdaTable || minKin << 372 if(1 < verboseLevel) {
238 << 373 G4cout << "### G4VEmProcess::BuildPhysicsTable() for "
239 G4EmTableUtil::BuildEmProcess(this, masterPr << 374 << GetProcessName()
240 nModels, verbo << 375 << " and particle " << num
241 isLocked, toBu << 376 << " buildLambdaTable= " << buildLambdaTable
>> 377 << " isTheMaster= " << isTheMaster
>> 378 << " " << masterProc
>> 379 << G4endl;
>> 380 }
>> 381
>> 382 if(particle == &part) {
>> 383
>> 384 // worker initialisation
>> 385 if(!isTheMaster) {
>> 386 theLambdaTable = masterProc->LambdaTable();
>> 387 theLambdaTablePrim = masterProc->LambdaTablePrim();
>> 388
>> 389 if(theLambdaTable) { FindLambdaMax(); }
>> 390
>> 391 // local initialisation of models
>> 392 G4bool printing = true;
>> 393 numberOfModels = modelManager->NumberOfModels();
>> 394 for(G4int i=0; i<numberOfModels; ++i) {
>> 395 G4VEmModel* mod = GetModelByIndex(i, printing);
>> 396 G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing);
>> 397 //G4cout << i << ". " << mod << " " << mod0 << " "
>> 398 // << particle->GetParticleName() << G4endl;
>> 399 mod->InitialiseLocal(particle, mod0);
>> 400 }
>> 401 // master thread
>> 402 } else {
>> 403 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
>> 404 BuildLambdaTable();
>> 405 }
>> 406 }
>> 407 }
>> 408
>> 409 // explicitly defined printout by particle name
>> 410 if(1 < verboseLevel ||
>> 411 (0 < verboseLevel && (num == "gamma" || num == "e-" ||
>> 412 num == "e+" || num == "mu+" ||
>> 413 num == "mu-" || num == "proton"||
>> 414 num == "pi+" || num == "pi-" ||
>> 415 num == "kaon+" || num == "kaon-" ||
>> 416 num == "alpha" || num == "anti_proton" ||
>> 417 num == "GenericIon"|| num == "alpha++" ||
>> 418 num == "alpha+" || num == "helium" ||
>> 419 num == "hydrogen")))
>> 420 {
>> 421 StreamInfo(G4cout, part);
>> 422 }
>> 423
>> 424 if(1 < verboseLevel) {
>> 425 G4cout << "### G4VEmProcess::BuildPhysicsTable() done for "
>> 426 << GetProcessName()
>> 427 << " and particle " << num
>> 428 << G4endl;
>> 429 }
242 } 430 }
243 431
244 //....oooOO0OOooo........oooOO0OOooo........oo 432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 433
246 void G4VEmProcess::BuildLambdaTable() 434 void G4VEmProcess::BuildLambdaTable()
247 { 435 {
>> 436 if(1 < verboseLevel) {
>> 437 G4cout << "G4EmProcess::BuildLambdaTable() for process "
>> 438 << GetProcessName() << " and particle "
>> 439 << particle->GetParticleName() << " " << this
>> 440 << G4endl;
>> 441 }
>> 442
>> 443 // Access to materials
>> 444 const G4ProductionCutsTable* theCoupleTable=
>> 445 G4ProductionCutsTable::GetProductionCutsTable();
>> 446 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 447
>> 448 G4LossTableBuilder* bld = lManager->GetTableBuilder();
>> 449
>> 450 G4PhysicsLogVector* aVector = nullptr;
>> 451 G4PhysicsLogVector* aVectorPrim = nullptr;
>> 452 G4PhysicsLogVector* bVectorPrim = nullptr;
>> 453
248 G4double scale = theParameters->MaxKinEnergy 454 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
249 G4int nbin = 455 G4int nbin =
250 theParameters->NumberOfBinsPerDecade()*G4l 456 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale));
>> 457 scale = G4Log(scale);
251 if(actBinning) { nbin = std::max(nbin, nLamb 458 if(actBinning) { nbin = std::max(nbin, nLambdaBins); }
252 scale = nbin/G4Log(scale); << 459 G4double emax1 = std::min(maxKinEnergy, minKinEnergyPrim);
253 << 460
254 G4LossTableBuilder* bld = lManager->GetTable << 461 for(size_t i=0; i<numOfCouples; ++i) {
255 G4EmTableUtil::BuildLambdaTable(this, partic << 462
256 bld, theLamb << 463 if (bld->GetFlag(i)) {
257 minKinEnergy << 464
258 maxKinEnergy << 465 // create physics vector and fill it
259 startFromNul << 466 const G4MaterialCutsCouple* couple =
>> 467 theCoupleTable->GetMaterialCutsCouple(i);
>> 468
>> 469 // build main table
>> 470 if(buildLambdaTable) {
>> 471 delete (*theLambdaTable)[i];
>> 472
>> 473 // if start from zero then change the scale
>> 474 G4double emin = minKinEnergy;
>> 475 G4bool startNull = false;
>> 476 if(startFromNull) {
>> 477 G4double e = MinPrimaryEnergy(particle,couple->GetMaterial());
>> 478 if(e >= emin) {
>> 479 emin = e;
>> 480 startNull = true;
>> 481 }
>> 482 }
>> 483 G4double emax = emax1;
>> 484 if(emax <= emin) { emax = 2*emin; }
>> 485 G4int bin = G4lrint(nbin*G4Log(emax/emin)/scale);
>> 486 if(bin < 3) { bin = 3; }
>> 487 aVector = new G4PhysicsLogVector(emin, emax, bin);
>> 488 aVector->SetSpline(splineFlag);
>> 489 modelManager->FillLambdaVector(aVector, couple, startNull);
>> 490 if(splineFlag) { aVector->FillSecondDerivatives(); }
>> 491 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
>> 492 }
>> 493 // build high energy table
>> 494 if(minKinEnergyPrim < maxKinEnergy) {
>> 495 delete (*theLambdaTablePrim)[i];
>> 496
>> 497 // start not from zero
>> 498 if(!bVectorPrim) {
>> 499 G4int bin = G4lrint(nbin*G4Log(maxKinEnergy/minKinEnergyPrim)/scale);
>> 500 if(bin < 3) { bin = 3; }
>> 501 aVectorPrim =
>> 502 new G4PhysicsLogVector(minKinEnergyPrim, maxKinEnergy, bin);
>> 503 bVectorPrim = aVectorPrim;
>> 504 } else {
>> 505 aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
>> 506 }
>> 507 // always use spline
>> 508 aVectorPrim->SetSpline(splineFlag);
>> 509 modelManager->FillLambdaVector(aVectorPrim, couple, false,
>> 510 fIsCrossSectionPrim);
>> 511 aVectorPrim->FillSecondDerivatives();
>> 512 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTablePrim, i,
>> 513 aVectorPrim);
>> 514 }
>> 515 }
>> 516 }
>> 517
>> 518 if(buildLambdaTable) { FindLambdaMax(); }
>> 519
>> 520 if(1 < verboseLevel) {
>> 521 G4cout << "Lambda table is built for "
>> 522 << particle->GetParticleName()
>> 523 << G4endl;
>> 524 }
260 } 525 }
261 526
262 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 528
264 void G4VEmProcess::StreamInfo(std::ostream& ou 529 void G4VEmProcess::StreamInfo(std::ostream& out,
265 const G4ParticleDefinition& 530 const G4ParticleDefinition& part, G4bool rst) const
266 { 531 {
267 G4String indent = (rst ? " " : ""); 532 G4String indent = (rst ? " " : "");
268 out << std::setprecision(6); 533 out << std::setprecision(6);
269 out << G4endl << indent << GetProcessName() 534 out << G4endl << indent << GetProcessName() << ": ";
270 if (!rst) { 535 if (!rst) {
271 out << " for " << part.GetParticleName(); 536 out << " for " << part.GetParticleName();
>> 537 if (integral) { out << ","; }
272 } 538 }
273 if(fXSType != fEmNoIntegral) { out << " XSt << 539 if(integral) { out << " integral:1 "; }
274 if(applyCuts) { out << " applyCuts:1 "; } 540 if(applyCuts) { out << " applyCuts:1 "; }
275 G4int subtype = GetProcessSubType(); << 541 out << " SubType=" << GetProcessSubType();
276 out << " SubType=" << subtype; << 542 if(biasFactor != 1.0) { out << " BiasingFactor= " << biasFactor; }
277 if (subtype == fAnnihilation) { <<
278 G4int mod = theParameters->PositronAtRestM <<
279 const G4String namp[2] = {"Simple", "Allis <<
280 out << " AtRestModel:" << namp[mod]; <<
281 } <<
282 if(biasFactor != 1.0) { out << " BiasingFac <<
283 out << " BuildTable=" << buildLambdaTable << 543 out << " BuildTable=" << buildLambdaTable << G4endl;
284 if(buildLambdaTable) { 544 if(buildLambdaTable) {
285 if(particle == &part) { 545 if(particle == &part) {
286 for(auto & v : *theLambdaTable) { << 546 size_t length = theLambdaTable->length();
287 if(nullptr != v) { << 547 for(size_t i=0; i<length; ++i) {
>> 548 G4PhysicsVector* v = (*theLambdaTable)[i];
>> 549 if(v) {
288 out << " Lambda table from "; 550 out << " Lambda table from ";
289 G4double emin = v->Energy(0); 551 G4double emin = v->Energy(0);
290 G4double emax = v->GetMaxEnergy(); 552 G4double emax = v->GetMaxEnergy();
291 G4int nbin = G4int(v->GetVectorLengt << 553 G4int nbin = v->GetVectorLength() - 1;
292 if(emin > minKinEnergy) { out << "th 554 if(emin > minKinEnergy) { out << "threshold "; }
293 else { out << G4BestUnit(emin,"Energ 555 else { out << G4BestUnit(emin,"Energy"); }
294 out << " to " 556 out << " to "
295 << G4BestUnit(emax,"Energy") 557 << G4BestUnit(emax,"Energy")
296 << ", " << G4lrint(nbin/std::log 558 << ", " << G4lrint(nbin/std::log10(emax/emin))
297 << " bins/decade, spline: " 559 << " bins/decade, spline: "
298 << splineFlag << G4endl; 560 << splineFlag << G4endl;
299 break; 561 break;
300 } 562 }
301 } 563 }
302 } else { 564 } else {
303 out << " Used Lambda table of " 565 out << " Used Lambda table of "
304 << particle->GetParticleName() << G4endl 566 << particle->GetParticleName() << G4endl;
305 } 567 }
306 } 568 }
307 if(minKinEnergyPrim < maxKinEnergy) { 569 if(minKinEnergyPrim < maxKinEnergy) {
308 if(particle == &part) { << 570 if(particle == &part) {
309 for(auto & v : *theLambdaTablePrim) { << 571 size_t length = theLambdaTablePrim->length();
310 if(nullptr != v) { << 572 for(size_t i=0; i<length; ++i) {
>> 573 G4PhysicsVector* v = (*theLambdaTablePrim)[i];
>> 574 if(v) {
311 out << " LambdaPrime table from 575 out << " LambdaPrime table from "
312 << G4BestUnit(v->Energy(0),"Ener 576 << G4BestUnit(v->Energy(0),"Energy")
313 << " to " 577 << " to "
314 << G4BestUnit(v->GetMaxEnergy(), 578 << G4BestUnit(v->GetMaxEnergy(),"Energy")
315 << " in " << v->GetVectorLength( 579 << " in " << v->GetVectorLength()-1
316 << " bins " << G4endl; 580 << " bins " << G4endl;
317 break; 581 break;
318 } 582 }
319 } 583 }
320 } else { 584 } else {
321 out << " Used LambdaPrime table of 585 out << " Used LambdaPrime table of "
322 << particle->GetParticleName() 586 << particle->GetParticleName() << G4endl;
323 } 587 }
324 } 588 }
325 StreamProcessInfo(out); 589 StreamProcessInfo(out);
326 modelManager->DumpModelList(out, verboseLeve 590 modelManager->DumpModelList(out, verboseLevel);
327 591
328 if(verboseLevel > 2 && buildLambdaTable) { 592 if(verboseLevel > 2 && buildLambdaTable) {
329 out << " LambdaTable address= " << th 593 out << " LambdaTable address= " << theLambdaTable << G4endl;
330 if(theLambdaTable && particle == &part) { 594 if(theLambdaTable && particle == &part) {
331 out << (*theLambdaTable) << G4endl; 595 out << (*theLambdaTable) << G4endl;
332 } 596 }
333 } 597 }
334 } 598 }
335 599
336 //....oooOO0OOooo........oooOO0OOooo........oo 600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337 601
338 void G4VEmProcess::StartTracking(G4Track* trac 602 void G4VEmProcess::StartTracking(G4Track* track)
339 { 603 {
340 // reset parameters for the new track 604 // reset parameters for the new track
341 currentParticle = track->GetParticleDefiniti 605 currentParticle = track->GetParticleDefinition();
342 theNumberOfInteractionLengthLeft = -1.0; 606 theNumberOfInteractionLengthLeft = -1.0;
343 mfpKinEnergy = DBL_MAX; 607 mfpKinEnergy = DBL_MAX;
344 preStepLambda = 0.0; <<
345 608
346 if(isIon) { massRatio = proton_mass_c2/curre 609 if(isIon) { massRatio = proton_mass_c2/currentParticle->GetPDGMass(); }
347 610
348 // forced biasing only for primary particles 611 // forced biasing only for primary particles
349 if(biasManager) { 612 if(biasManager) {
350 if(0 == track->GetParentID()) { 613 if(0 == track->GetParentID()) {
351 // primary particle 614 // primary particle
352 biasFlag = true; 615 biasFlag = true;
353 biasManager->ResetForcedInteraction(); 616 biasManager->ResetForcedInteraction();
354 } 617 }
355 } 618 }
356 } 619 }
357 620
358 //....oooOO0OOooo........oooOO0OOooo........oo 621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 622
360 G4double G4VEmProcess::PostStepGetPhysicalInte 623 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& tr 624 const G4Track& track,
362 G4double previo 625 G4double previousStepSize,
363 G4ForceCondition* 626 G4ForceCondition* condition)
364 { 627 {
365 *condition = NotForced; 628 *condition = NotForced;
366 G4double x = DBL_MAX; 629 G4double x = DBL_MAX;
367 630
368 DefineMaterial(track.GetMaterialCutsCouple() 631 DefineMaterial(track.GetMaterialCutsCouple());
369 preStepKinEnergy = track.GetKineticEnergy(); << 632 preStepKinEnergy = track.GetKineticEnergy();
370 const G4double scaledEnergy = preStepKinEner << 633 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
>> 634 G4double scaledEnergy = preStepKinEnergy*massRatio;
371 SelectModel(scaledEnergy, currentCoupleIndex 635 SelectModel(scaledEnergy, currentCoupleIndex);
372 /* << 636
373 G4cout << "PostStepGetPhysicalInteractionLen <<
374 << " couple: " << currentCouple << G <<
375 */ <<
376 if(!currentModel->IsActive(scaledEnergy)) { 637 if(!currentModel->IsActive(scaledEnergy)) {
377 theNumberOfInteractionLengthLeft = -1.0; 638 theNumberOfInteractionLengthLeft = -1.0;
378 currentInteractionLength = DBL_MAX; 639 currentInteractionLength = DBL_MAX;
379 mfpKinEnergy = DBL_MAX; <<
380 preStepLambda = 0.0; <<
381 return x; 640 return x;
382 } 641 }
383 642
384 // forced biasing only for primary particles 643 // forced biasing only for primary particles
385 if(biasManager) { 644 if(biasManager) {
386 if(0 == track.GetParentID()) { 645 if(0 == track.GetParentID()) {
387 if(biasFlag && 646 if(biasFlag &&
388 biasManager->ForcedInteractionRegion( << 647 biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
389 return biasManager->GetStepLimit((G4in << 648 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
390 } 649 }
391 } 650 }
392 } 651 }
393 652
394 // compute mean free path 653 // compute mean free path
>> 654 if(preStepKinEnergy < mfpKinEnergy) {
>> 655 if (integral) {
>> 656 ComputeIntegralLambda(preStepKinEnergy, preStepLogKinEnergy);
>> 657 } else {
>> 658 preStepLambda = GetCurrentLambda(preStepKinEnergy, preStepLogKinEnergy);
>> 659 }
395 660
396 ComputeIntegralLambda(preStepKinEnergy, trac << 661 // zero cross section
397 << 662 if(preStepLambda <= 0.0) {
398 // zero cross section << 663 theNumberOfInteractionLengthLeft = -1.0;
399 if(preStepLambda <= 0.0) { << 664 currentInteractionLength = DBL_MAX;
400 theNumberOfInteractionLengthLeft = -1.0; << 665 }
401 currentInteractionLength = DBL_MAX; << 666 }
402 667
403 } else { << 668 // non-zero cross section
>> 669 if(preStepLambda > 0.0) {
404 670
405 // non-zero cross section <<
406 if (theNumberOfInteractionLengthLeft < 0.0 671 if (theNumberOfInteractionLengthLeft < 0.0) {
407 672
408 // beggining of tracking (or just after 673 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Lo << 674 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
410 theInitialNumberOfInteractionLength = th 675 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
411 676
412 } else { << 677 } else if(currentInteractionLength < DBL_MAX) {
413 678
414 theNumberOfInteractionLengthLeft -= 679 theNumberOfInteractionLengthLeft -=
415 previousStepSize/currentInteractionLen 680 previousStepSize/currentInteractionLength;
416 theNumberOfInteractionLengthLeft = 681 theNumberOfInteractionLengthLeft =
417 std::max(theNumberOfInteractionLengthL 682 std::max(theNumberOfInteractionLengthLeft, 0.0);
418 } 683 }
419 684
420 // new mean free path and step limit for t 685 // new mean free path and step limit for the next step
421 currentInteractionLength = 1.0/preStepLamb 686 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * cur 687 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 } 688 }
424 return x; 689 return x;
425 } 690 }
426 691
427 //....oooOO0OOooo........oooOO0OOooo........oo 692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428 693
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 694 void G4VEmProcess::ComputeIntegralLambda(G4double e, G4double loge)
430 { 695 {
431 if (fXSType == fEmNoIntegral) { << 696 // condition to skip recomputation of cross section
432 preStepLambda = GetCurrentLambda(e, LogEki << 697 const G4double epeak = theEnergyOfCrossSectionMax[currentCoupleIndex];
433 << 698 if(e <= epeak && e/lambdaFactor >= mfpKinEnergy) { return; }
434 } else if (fXSType == fEmIncreasing) { << 699
435 if(e*invLambdaFactor < mfpKinEnergy) { << 700 // recomputation is needed
436 preStepLambda = GetCurrentLambda(e, LogE << 701 if (e <= epeak) {
437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 702 preStepLambda = GetCurrentLambda(e, loge);
438 } << 703 mfpKinEnergy = e;
439 <<
440 } else if(fXSType == fEmDecreasing) { <<
441 if(e < mfpKinEnergy) { <<
442 const G4double e1 = e*lambdaFactor; <<
443 preStepLambda = GetCurrentLambda(e1); <<
444 mfpKinEnergy = e1; <<
445 } <<
446 <<
447 } else if(fXSType == fEmOnePeak) { <<
448 const G4double epeak = (*theEnergyOfCrossS <<
449 if(e <= epeak) { <<
450 if(e*invLambdaFactor < mfpKinEnergy) { <<
451 preStepLambda = GetCurrentLambda(e, Lo <<
452 mfpKinEnergy = (preStepLambda > 0.0) ? <<
453 } <<
454 } else if(e < mfpKinEnergy) { <<
455 const G4double e1 = std::max(epeak, e*la <<
456 preStepLambda = GetCurrentLambda(e1); <<
457 mfpKinEnergy = e1; <<
458 } <<
459 } else { 704 } else {
460 preStepLambda = GetCurrentLambda(e, LogEki << 705 const G4double e1 = e*lambdaFactor;
>> 706 if (e1 > epeak) {
>> 707 preStepLambda = GetCurrentLambda(e, loge);
>> 708 mfpKinEnergy = e;
>> 709 const G4double preStepLambda1 = GetCurrentLambda(e1,loge+logLambdaFactor);
>> 710 if (preStepLambda1 > preStepLambda) {
>> 711 mfpKinEnergy = e1;
>> 712 preStepLambda = preStepLambda1;
>> 713 }
>> 714 } else {
>> 715 preStepLambda = fFactor*theCrossSectionMax[currentCoupleIndex];
>> 716 mfpKinEnergy = epeak;
>> 717 }
461 } 718 }
462 } 719 }
463 720
464 //....oooOO0OOooo........oooOO0OOooo........oo 721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 722
466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 723 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 724 const G4Step& step)
468 { 725 {
469 // clear number of interaction lengths in an << 726 // In all cases clear number of interaction lengths
470 theNumberOfInteractionLengthLeft = -1.0; 727 theNumberOfInteractionLengthLeft = -1.0;
471 mfpKinEnergy = DBL_MAX; << 728 mfpKinEnergy = DBL_MAX;
472 729
473 fParticleChange.InitializeForPostStep(track) 730 fParticleChange.InitializeForPostStep(track);
474 731
475 // Do not make anything if particle is stopp 732 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt! 733 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) 734 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478 735
479 const G4double finalT = track.GetKineticEner << 736 const G4double finalT = track.GetKineticEnergy();
>> 737 const G4double logFinalT = track.GetDynamicParticle()->GetLogKineticEnergy();
480 738
481 // forced process - should happen only once 739 // forced process - should happen only once per track
482 if(biasFlag) { 740 if(biasFlag) {
483 if(biasManager->ForcedInteractionRegion((G << 741 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
484 biasFlag = false; 742 biasFlag = false;
485 } 743 }
486 } 744 }
487 745
488 // check active and select model <<
489 const G4double scaledEnergy = finalT*massRat <<
490 SelectModel(scaledEnergy, currentCoupleIndex <<
491 if(!currentModel->IsActive(scaledEnergy)) { <<
492 <<
493 // Integral approach 746 // Integral approach
494 if (fXSType != fEmNoIntegral) { << 747 if (integral) {
495 const G4double logFinalT = << 748 G4double lx = GetLambda(finalT, currentCouple, logFinalT);
496 track.GetDynamicParticle()->GetLogKineti << 749 if(preStepLambda<lx && 1 < verboseLevel) {
497 const G4double lx = std::max(GetCurrentLam <<
498 #ifdef G4VERBOSE <<
499 if(preStepLambda < lx && 1 < verboseLevel) <<
500 G4cout << "WARNING: for " << currentPart 750 G4cout << "WARNING: for " << currentParticle->GetParticleName()
501 << " and " << GetProcessName() << << 751 << " and " << GetProcessName()
502 << " preLambda= " << preStepLambd << 752 << " E(MeV)= " << finalT/MeV
503 << " < " << lx << " (postLambda) << 753 << " preLambda= " << preStepLambda << " < "
504 } << 754 << lx << " (postLambda) "
505 #endif << 755 << G4endl;
506 // if false interaction then use new cross << 756 }
507 // if both values are zero - no interactio << 757
508 if(preStepLambda*G4UniformRand() >= lx) { << 758 if(preStepLambda*G4UniformRand() > lx) {
>> 759 ClearNumberOfInteractionLengthLeft();
509 return &fParticleChange; 760 return &fParticleChange;
510 } 761 }
511 } 762 }
512 763
>> 764 G4double scaledEnergy = finalT*massRatio;
>> 765 SelectModel(scaledEnergy, currentCoupleIndex);
>> 766 if(!currentModel->IsActive(scaledEnergy)) { return &fParticleChange; }
>> 767
513 // define new weight for primary and seconda 768 // define new weight for primary and secondaries
514 G4double weight = fParticleChange.GetParentW 769 G4double weight = fParticleChange.GetParentWeight();
515 if(weightFlag) { 770 if(weightFlag) {
516 weight /= biasFactor; 771 weight /= biasFactor;
517 fParticleChange.ProposeWeight(weight); 772 fParticleChange.ProposeWeight(weight);
518 } 773 }
519 774
520 #ifdef G4VERBOSE << 775
521 if(1 < verboseLevel) { 776 if(1 < verboseLevel) {
522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 777 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV 778 << finalT/MeV
524 << " MeV; model= (" << currentModel 779 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnerg 780 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl; 781 << G4endl;
527 } 782 }
528 #endif << 783
529 784
530 // sample secondaries 785 // sample secondaries
531 secParticles.clear(); 786 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticle 787 currentModel->SampleSecondaries(&secParticles,
533 currentCoupl 788 currentCouple,
534 track.GetDyn 789 track.GetDynamicParticle(),
535 (*theCuts)[c 790 (*theCuts)[currentCoupleIndex]);
536 791
537 G4int num0 = (G4int)secParticles.size(); << 792 G4int num0 = secParticles.size();
538 793
539 // splitting or Russian roulette 794 // splitting or Russian roulette
540 if(biasManager) { 795 if(biasManager) {
541 if(biasManager->SecondaryBiasingRegion((G4 << 796 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
542 G4double eloss = 0.0; 797 G4double eloss = 0.0;
543 weight *= biasManager->ApplySecondaryBia 798 weight *= biasManager->ApplySecondaryBiasing(
544 secParticles, track, currentModel, &fP 799 secParticles, track, currentModel, &fParticleChange, eloss,
545 (G4int)currentCoupleIndex, (*theCuts)[ << 800 currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546 step.GetPostStepPoint()->GetSafety()); 801 step.GetPostStepPoint()->GetSafety());
547 if(eloss > 0.0) { 802 if(eloss > 0.0) {
548 eloss += fParticleChange.GetLocalEnerg 803 eloss += fParticleChange.GetLocalEnergyDeposit();
549 fParticleChange.ProposeLocalEnergyDepo 804 fParticleChange.ProposeLocalEnergyDeposit(eloss);
550 } 805 }
551 } 806 }
552 } 807 }
553 808
554 // save secondaries 809 // save secondaries
555 G4int num = (G4int)secParticles.size(); << 810 G4int num = secParticles.size();
556 if(num > 0) { 811 if(num > 0) {
557 812
558 fParticleChange.SetNumberOfSecondaries(num 813 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEn 814 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime(); 815 G4double time = track.GetGlobalTime();
561 <<
562 G4int n1(0), n2(0); <<
563 if(num0 > mainSecondaries) { <<
564 currentModel->FillNumberOfSecondaries(n1 <<
565 } <<
566 816
567 for (G4int i=0; i<num; ++i) { 817 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i]; << 818 if (secParticles[i]) {
569 if (nullptr != dp) { << 819 G4DynamicParticle* dp = secParticles[i];
570 const G4ParticleDefinition* p = dp->Ge 820 const G4ParticleDefinition* p = dp->GetParticleDefinition();
571 G4double e = dp->GetKineticEnergy(); 821 G4double e = dp->GetKineticEnergy();
572 G4bool good = true; 822 G4bool good = true;
573 if(applyCuts) { 823 if(applyCuts) {
574 if (p == theGamma) { 824 if (p == theGamma) {
575 if (e < (*theCutsGamma)[currentCou 825 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576 826
577 } else if (p == theElectron) { 827 } else if (p == theElectron) {
578 if (e < (*theCutsElectron)[current 828 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579 829
580 } else if (p == thePositron) { 830 } else if (p == thePositron) {
581 if (electron_mass_c2 < (*theCutsGa 831 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582 e < (*theCutsPositron)[current 832 e < (*theCutsPositron)[currentCoupleIndex]) {
583 good = false; 833 good = false;
584 e += 2.0*electron_mass_c2; 834 e += 2.0*electron_mass_c2;
585 } 835 }
586 } 836 }
587 // added secondary if it is good 837 // added secondary if it is good
588 } 838 }
589 if (good) { 839 if (good) {
590 G4Track* t = new G4Track(dp, time, t 840 G4Track* t = new G4Track(dp, time, track.GetPosition());
591 t->SetTouchableHandle(track.GetTouch 841 t->SetTouchableHandle(track.GetTouchableHandle());
592 if (biasManager) { 842 if (biasManager) {
593 t->SetWeight(weight * biasManager- 843 t->SetWeight(weight * biasManager->GetWeight(i));
594 } else { 844 } else {
595 t->SetWeight(weight); 845 t->SetWeight(weight);
596 } 846 }
597 pParticleChange->AddSecondary(t); 847 pParticleChange->AddSecondary(t);
598 848
599 // define type of secondary 849 // define type of secondary
600 if(i < mainSecondaries) { << 850 if(i < mainSecondaries) { t->SetCreatorModelIndex(secID); }
601 t->SetCreatorModelID(secID); << 851 else if(i < num0) {
602 if(GetProcessSubType() == fCompton << 852 if(p == theGamma) {
603 t->SetCreatorModelID(_ComptonGam << 853 t->SetCreatorModelIndex(fluoID);
604 } <<
605 } else if(i < mainSecondaries + n1) <<
606 t->SetCreatorModelID(tripletID); <<
607 } else if(i < mainSecondaries + n1 + <<
608 t->SetCreatorModelID(_IonRecoil); <<
609 } else { <<
610 if(i < num0) { <<
611 if(p == theGamma) { <<
612 t->SetCreatorModelID(fluoID); <<
613 } else { <<
614 t->SetCreatorModelID(augerID); <<
615 } <<
616 } else { 854 } else {
617 t->SetCreatorModelID(biasID); << 855 t->SetCreatorModelIndex(augerID);
618 } 856 }
>> 857 } else {
>> 858 t->SetCreatorModelIndex(biasID);
619 } 859 }
620 /* 860 /*
621 G4cout << "Secondary(post step) has 861 G4cout << "Secondary(post step) has weight " << t->GetWeight()
622 << ", Ekin= " << t->GetKineti 862 << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV "
623 << GetProcessName() << " fluo 863 << GetProcessName() << " fluoID= " << fluoID
624 << " augerID= " << augerID << 864 << " augerID= " << augerID <<G4endl;
625 */ 865 */
626 } else { 866 } else {
627 delete dp; 867 delete dp;
628 edep += e; 868 edep += e;
629 } 869 }
630 } 870 }
631 } 871 }
632 fParticleChange.ProposeLocalEnergyDeposit( 872 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 } 873 }
634 874
635 if(0.0 == fParticleChange.GetProposedKinetic 875 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636 fAlive == fParticleChange.GetTrackStatus( 876 fAlive == fParticleChange.GetTrackStatus()) {
637 if(particle->GetProcessManager()->GetAtRes 877 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638 { fParticleChange.ProposeTrackStatus( 878 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639 else { fParticleChange.ProposeTrackStatus( 879 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640 } 880 }
641 881
642 return &fParticleChange; 882 return &fParticleChange;
643 } 883 }
644 884
645 //....oooOO0OOooo........oooOO0OOooo........oo 885 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 886
647 G4bool G4VEmProcess::StorePhysicsTable(const G 887 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G 888 const G4String& directory,
649 G4bool 889 G4bool ascii)
650 { 890 {
651 if(!isTheMaster || part != particle) { retur << 891 G4bool yes = true;
652 if(G4EmTableUtil::StoreTable(this, part, the << 892 if(!isTheMaster) { return yes; }
653 directory, "Lambda", << 893
654 verboseLevel, a << 894 if ( theLambdaTable && part == particle) {
655 G4EmTableUtil::StoreTable(this, part, the << 895 const G4String& nam =
656 directory, "LambdaPrim", << 896 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
657 verboseLevel, a << 897 yes = theLambdaTable->StorePhysicsTable(nam,ascii);
658 return true; << 898
>> 899 if ( yes ) {
>> 900 if(0 < verboseLevel) G4cout << "Stored: " << nam << G4endl;
>> 901 } else {
>> 902 G4cout << "Fail to store Physics Table for "
>> 903 << particle->GetParticleName()
>> 904 << " and process " << GetProcessName()
>> 905 << " in the directory <" << directory
>> 906 << "> " << G4endl;
>> 907 }
659 } 908 }
660 return false; << 909 if ( theLambdaTablePrim && part == particle) {
>> 910 const G4String& name =
>> 911 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
>> 912 yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
>> 913
>> 914 if ( yes ) {
>> 915 if(0 < verboseLevel) {
>> 916 G4cout << "Physics table prim is stored for "
>> 917 << particle->GetParticleName()
>> 918 << " and process " << GetProcessName()
>> 919 << " in the directory <" << directory
>> 920 << "> " << G4endl;
>> 921 }
>> 922 } else {
>> 923 G4cout << "Fail to store Physics Table Prim for "
>> 924 << particle->GetParticleName()
>> 925 << " and process " << GetProcessName()
>> 926 << " in the directory <" << directory
>> 927 << "> " << G4endl;
>> 928 }
>> 929 }
>> 930 return yes;
661 } 931 }
662 932
663 //....oooOO0OOooo........oooOO0OOooo........oo 933 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664 934
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 935 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 cons << 936 const G4String& directory,
667 G4bo 937 G4bool ascii)
668 { 938 {
669 if(!isTheMaster || part != particle) { retur << 939 if(1 < verboseLevel) {
>> 940 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
>> 941 << part->GetParticleName() << " and process "
>> 942 << GetProcessName() << G4endl;
>> 943 }
670 G4bool yes = true; 944 G4bool yes = true;
>> 945
>> 946 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy)
>> 947 || particle != part) { return yes; }
>> 948
>> 949 const G4String particleName = part->GetParticleName();
>> 950
671 if(buildLambdaTable) { 951 if(buildLambdaTable) {
672 yes = G4EmTableUtil::RetrieveTable(this, p << 952 const G4String& filename =
673 "Lambda << 953 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
674 ascii, << 954 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
675 } << 955 filename,ascii);
676 if(yes && minKinEnergyPrim < maxKinEnergy) { << 956 if ( yes ) {
677 yes = G4EmTableUtil::RetrieveTable(this, p << 957 if (0 < verboseLevel) {
678 "Lambda << 958 G4cout << "Lambda table for " << particleName
679 ascii, << 959 << " is Retrieved from <"
>> 960 << filename << ">"
>> 961 << G4endl;
>> 962 }
>> 963 if(theParameters->Spline()) {
>> 964 size_t n = theLambdaTable->length();
>> 965 for(size_t i=0; i<n; ++i) {
>> 966 if((* theLambdaTable)[i]) {
>> 967 (* theLambdaTable)[i]->SetSpline(true);
>> 968 }
>> 969 }
>> 970 }
>> 971 } else {
>> 972 if (1 < verboseLevel) {
>> 973 G4cout << "Lambda table for " << particleName << " in file <"
>> 974 << filename << "> is not exist"
>> 975 << G4endl;
>> 976 }
>> 977 }
>> 978 }
>> 979 if(minKinEnergyPrim < maxKinEnergy) {
>> 980 const G4String& filename =
>> 981 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
>> 982 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTablePrim,
>> 983 filename,ascii);
>> 984 if ( yes ) {
>> 985 if (0 < verboseLevel) {
>> 986 G4cout << "Lambda table prim for " << particleName
>> 987 << " is Retrieved from <"
>> 988 << filename << ">"
>> 989 << G4endl;
>> 990 }
>> 991 if(theParameters->Spline()) {
>> 992 size_t n = theLambdaTablePrim->length();
>> 993 for(size_t i=0; i<n; ++i) {
>> 994 if((* theLambdaTablePrim)[i]) {
>> 995 (* theLambdaTablePrim)[i]->SetSpline(true);
>> 996 }
>> 997 }
>> 998 }
>> 999 } else {
>> 1000 if (1 < verboseLevel) {
>> 1001 G4cout << "Lambda table prim for " << particleName << " in file <"
>> 1002 << filename << "> is not exist"
>> 1003 << G4endl;
>> 1004 }
>> 1005 }
680 } 1006 }
681 return yes; 1007 return yes;
682 } 1008 }
683 1009
684 //....oooOO0OOooo........oooOO0OOooo........oo 1010 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685 1011
686 G4double G4VEmProcess::GetCrossSection(G4doubl << 1012 G4double
687 const G << 1013 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
>> 1014 const G4MaterialCutsCouple* couple,
>> 1015 G4double logKinEnergy)
688 { 1016 {
689 CurrentSetup(couple, kinEnergy); << 1017 // Cross section per atom is calculated
690 return GetCurrentLambda(kinEnergy, G4Log(kin << 1018 DefineMaterial(couple);
>> 1019 G4double cross = 0.0;
>> 1020 if(buildLambdaTable) {
>> 1021 cross = GetCurrentLambda(kineticEnergy,
>> 1022 (logKinEnergy < DBL_MAX) ? logKinEnergy : G4Log(kineticEnergy));
>> 1023 } else {
>> 1024 SelectModel(kineticEnergy, currentCoupleIndex);
>> 1025 if(currentModel) {
>> 1026 cross = fFactor*currentModel->CrossSectionPerVolume(currentMaterial,
>> 1027 currentParticle,
>> 1028 kineticEnergy);
>> 1029 }
>> 1030 }
>> 1031 return std::max(cross, 0.0);
691 } 1032 }
692 1033
693 //....oooOO0OOooo........oooOO0OOooo........oo 1034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694 1035
695 G4double G4VEmProcess::GetMeanFreePath(const G 1036 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696 G4doubl 1037 G4double,
697 G4Force 1038 G4ForceCondition* condition)
698 { 1039 {
699 *condition = NotForced; 1040 *condition = NotForced;
700 return G4VEmProcess::MeanFreePath(track); 1041 return G4VEmProcess::MeanFreePath(track);
701 } 1042 }
702 1043
703 //....oooOO0OOooo........oooOO0OOooo........oo 1044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 1045
>> 1046 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
>> 1047 {
>> 1048 const G4double kinEnergy = track.GetKineticEnergy();
>> 1049 CurrentSetup(track.GetMaterialCutsCouple(), kinEnergy);
>> 1050 const G4double xs = GetCurrentLambda(kinEnergy,
>> 1051 track.GetDynamicParticle()->GetLogKineticEnergy());
>> 1052 return (0.0 < xs) ? 1.0/xs : DBL_MAX;
>> 1053 }
>> 1054
>> 1055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 1056
705 G4double 1057 G4double
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou 1058 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kinEnergy,
707 G4dou 1059 G4double Z, G4double A, G4double cut)
708 { 1060 {
709 SelectModel(kinEnergy, currentCoupleIndex); 1061 SelectModel(kinEnergy, currentCoupleIndex);
710 return (currentModel) ? 1062 return (currentModel) ?
711 currentModel->ComputeCrossSectionPerAtom(c 1063 currentModel->ComputeCrossSectionPerAtom(currentParticle, kinEnergy,
712 Z 1064 Z, A, cut) : 0.0;
713 } 1065 }
714 1066
715 //....oooOO0OOooo........oooOO0OOooo........oo 1067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 1068
717 G4PhysicsVector* << 1069 void G4VEmProcess::FindLambdaMax()
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate <<
719 { 1070 {
720 DefineMaterial(couple); << 1071 if(1 < verboseLevel) {
721 G4PhysicsVector* newv = new G4PhysicsLogVect << 1072 G4cout << "### G4VEmProcess::FindLambdaMax: "
722 << 1073 << particle->GetParticleName()
723 return newv; << 1074 << " and process " << GetProcessName() << " " << G4endl;
724 } << 1075 }
725 << 1076 size_t n = theLambdaTable->length();
726 //....oooOO0OOooo........oooOO0OOooo........oo << 1077 G4PhysicsVector* pv;
727 << 1078 G4double e, ss, emax, smax;
728 const G4Element* G4VEmProcess::GetCurrentEleme << 1079
729 { << 1080 size_t i;
730 return (nullptr != currentModel) ? << 1081
731 currentModel->GetCurrentElement(currentMat << 1082 // first loop on existing vectors
>> 1083 for (i=0; i<n; ++i) {
>> 1084 pv = (*theLambdaTable)[i];
>> 1085 if(pv) {
>> 1086 size_t nb = pv->GetVectorLength();
>> 1087 emax = DBL_MAX;
>> 1088 smax = 0.0;
>> 1089 if(nb > 0) {
>> 1090 for (size_t j=0; j<nb; ++j) {
>> 1091 e = pv->Energy(j);
>> 1092 ss = (*pv)(j);
>> 1093 if(ss > smax) {
>> 1094 smax = ss;
>> 1095 emax = e;
>> 1096 }
>> 1097 }
>> 1098 }
>> 1099 theEnergyOfCrossSectionMax[i] = emax;
>> 1100 theCrossSectionMax[i] = smax;
>> 1101 if(1 < verboseLevel) {
>> 1102 G4cout << "For " << particle->GetParticleName()
>> 1103 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
>> 1104 << " lambda= " << smax << G4endl;
>> 1105 }
>> 1106 }
>> 1107 }
>> 1108 // second loop using base materials
>> 1109 for (i=0; i<n; ++i) {
>> 1110 pv = (*theLambdaTable)[i];
>> 1111 if(!pv){
>> 1112 G4int j = (*theDensityIdx)[i];
>> 1113 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
>> 1114 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
>> 1115 }
>> 1116 }
732 } 1117 }
733 1118
734 //....oooOO0OOooo........oooOO0OOooo........oo 1119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735 1120
736 const G4Element* G4VEmProcess::GetTargetElemen << 1121 G4PhysicsVector*
>> 1122 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
737 { 1123 {
738 return (nullptr != currentModel) ? << 1124 G4PhysicsVector* v =
739 currentModel->GetCurrentElement(currentMat << 1125 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
>> 1126 v->SetSpline(theParameters->Spline());
>> 1127 return v;
740 } 1128 }
741 1129
742 //....oooOO0OOooo........oooOO0OOooo........oo 1130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743 1131
744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 1132 const G4Element* G4VEmProcess::GetCurrentElement() const
745 { 1133 {
746 return (nullptr != currentModel) ? << 1134 const G4Element* elm =
747 currentModel->GetCurrentIsotope(GetCurrent << 1135 (currentModel) ? currentModel->GetCurrentElement() : nullptr;
>> 1136 return elm;
748 } 1137 }
749 1138
750 //....oooOO0OOooo........oooOO0OOooo........oo 1139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 1140
752 void G4VEmProcess::SetCrossSectionBiasingFacto 1141 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 { 1142 {
754 if(f > 0.0) { 1143 if(f > 0.0) {
755 biasFactor = f; 1144 biasFactor = f;
756 weightFlag = flag; 1145 weightFlag = flag;
757 if(1 < verboseLevel) { 1146 if(1 < verboseLevel) {
758 G4cout << "### SetCrossSectionBiasingFac 1147 G4cout << "### SetCrossSectionBiasingFactor: for "
759 << particle->GetParticleName() 1148 << particle->GetParticleName()
760 << " and process " << GetProcessN 1149 << " and process " << GetProcessName()
761 << " biasFactor= " << f << " weig 1150 << " biasFactor= " << f << " weightFlag= " << flag
762 << G4endl; 1151 << G4endl;
763 } 1152 }
764 } 1153 }
765 } 1154 }
766 1155
767 //....oooOO0OOooo........oooOO0OOooo........oo 1156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768 1157
769 void 1158 void
770 G4VEmProcess::ActivateForcedInteraction(G4doub 1159 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771 G4bool 1160 G4bool flag)
772 { 1161 {
773 if(nullptr == biasManager) { biasManager = n << 1162 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
774 if(1 < verboseLevel) { 1163 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: 1164 G4cout << "### ActivateForcedInteraction: for "
776 << particle->GetParticleName() 1165 << particle->GetParticleName()
777 << " and process " << GetProcessNam 1166 << " and process " << GetProcessName()
778 << " length(mm)= " << length/mm 1167 << " length(mm)= " << length/mm
779 << " in G4Region <" << r 1168 << " in G4Region <" << r
780 << "> weightFlag= " << flag 1169 << "> weightFlag= " << flag
781 << G4endl; 1170 << G4endl;
782 } 1171 }
783 weightFlag = flag; 1172 weightFlag = flag;
784 biasManager->ActivateForcedInteraction(lengt 1173 biasManager->ActivateForcedInteraction(length, r);
785 } 1174 }
786 1175
787 //....oooOO0OOooo........oooOO0OOooo........oo 1176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788 1177
789 void 1178 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G 1179 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791 G4double factor, 1180 G4double factor,
792 G4double energyLimit) 1181 G4double energyLimit)
793 { 1182 {
794 if (0.0 <= factor) { 1183 if (0.0 <= factor) {
795 1184
796 // Range cut can be applied only for e- 1185 // Range cut can be applied only for e-
797 if(0.0 == factor && secondaryParticle != G 1186 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798 { return; } 1187 { return; }
799 1188
800 if(!biasManager) { biasManager = new G4EmB 1189 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801 biasManager->ActivateSecondaryBiasing(regi 1190 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802 if(1 < verboseLevel) { 1191 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: 1192 G4cout << "### ActivateSecondaryBiasing: for "
804 << " process " << GetProcessName() 1193 << " process " << GetProcessName()
805 << " factor= " << factor 1194 << " factor= " << factor
806 << " in G4Region <" << region 1195 << " in G4Region <" << region
807 << "> energyLimit(MeV)= " << energyLimi 1196 << "> energyLimit(MeV)= " << energyLimit/MeV
808 << G4endl; 1197 << G4endl;
809 } 1198 }
810 } 1199 }
811 } 1200 }
812 1201
813 //....oooOO0OOooo........oooOO0OOooo........oo 1202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 1203
815 void G4VEmProcess::SetLambdaBinning(G4int n) 1204 void G4VEmProcess::SetLambdaBinning(G4int n)
816 { 1205 {
817 if(5 < n && n < 10000000) { 1206 if(5 < n && n < 10000000) {
818 nLambdaBins = n; 1207 nLambdaBins = n;
819 actBinning = true; 1208 actBinning = true;
820 } else { 1209 } else {
821 G4double e = (G4double)n; 1210 G4double e = (G4double)n;
822 PrintWarning("SetLambdaBinning", e); 1211 PrintWarning("SetLambdaBinning", e);
823 } 1212 }
824 } 1213 }
825 1214
826 //....oooOO0OOooo........oooOO0OOooo........oo 1215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827 1216
828 void G4VEmProcess::SetMinKinEnergy(G4double e) 1217 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 { 1218 {
830 if(1.e-3*eV < e && e < maxKinEnergy) { 1219 if(1.e-3*eV < e && e < maxKinEnergy) {
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma 1220 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e)
832 /G4Log(maxKinEnergy/ 1221 /G4Log(maxKinEnergy/minKinEnergy));
833 minKinEnergy = e; 1222 minKinEnergy = e;
834 actMinKinEnergy = true; 1223 actMinKinEnergy = true;
835 } else { PrintWarning("SetMinKinEnergy", e); 1224 } else { PrintWarning("SetMinKinEnergy", e); }
836 } 1225 }
837 1226
838 //....oooOO0OOooo........oooOO0OOooo........oo 1227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839 1228
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) 1229 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 { 1230 {
842 if(minKinEnergy < e && e < 1.e+6*TeV) { 1231 if(minKinEnergy < e && e < 1.e+6*TeV) {
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ 1232 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy)
844 /G4Log(maxKinEnergy/ 1233 /G4Log(maxKinEnergy/minKinEnergy));
845 maxKinEnergy = e; 1234 maxKinEnergy = e;
846 actMaxKinEnergy = true; 1235 actMaxKinEnergy = true;
847 } else { PrintWarning("SetMaxKinEnergy", e); 1236 } else { PrintWarning("SetMaxKinEnergy", e); }
848 } 1237 }
849 1238
850 //....oooOO0OOooo........oooOO0OOooo........oo 1239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851 1240
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl 1241 void G4VEmProcess::SetMinKinEnergyPrim(G4double e)
853 { 1242 {
854 if(theParameters->MinKinEnergy() <= e && 1243 if(theParameters->MinKinEnergy() <= e &&
855 e <= theParameters->MaxKinEnergy()) { min 1244 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; }
856 else { PrintWarning("SetMinKinEnergyPrim", e 1245 else { PrintWarning("SetMinKinEnergyPrim", e); }
857 } 1246 }
858 1247
859 //....oooOO0OOooo........oooOO0OOooo........oo 1248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
860 1249
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const 1250 G4VEmProcess* G4VEmProcess::GetEmProcess(const G4String& nam)
862 { 1251 {
863 return (nam == GetProcessName()) ? this : nu 1252 return (nam == GetProcessName()) ? this : nullptr;
864 } 1253 }
865 1254
866 //....oooOO0OOooo........oooOO0OOooo........oo 1255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
867 1256
868 G4double G4VEmProcess::PolarAngleLimit() const << 1257 G4double
>> 1258 G4VEmProcess::GetLambda(G4double kinEnergy, const G4MaterialCutsCouple* couple)
869 { 1259 {
870 return theParameters->MscThetaLimit(); << 1260 CurrentSetup(couple, kinEnergy);
>> 1261 return GetCurrentLambda(kinEnergy, G4Log(kinEnergy));
871 } 1262 }
872 1263
873 //....oooOO0OOooo........oooOO0OOooo........oo 1264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
874 1265
875 void G4VEmProcess::PrintWarning(G4String tit, 1266 void G4VEmProcess::PrintWarning(G4String tit, G4double val)
876 { 1267 {
877 G4String ss = "G4VEmProcess::" + tit; 1268 G4String ss = "G4VEmProcess::" + tit;
878 G4ExceptionDescription ed; 1269 G4ExceptionDescription ed;
879 ed << "Parameter is out of range: " << val 1270 ed << "Parameter is out of range: " << val
880 << " it will have no effect!\n" << " Pro 1271 << " it will have no effect!\n" << " Process "
881 << GetProcessName() << " nbins= " << the 1272 << GetProcessName() << " nbins= " << theParameters->NumberOfBins()
882 << " Emin(keV)= " << theParameters->MinKi 1273 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV
883 << " Emax(GeV)= " << theParameters->MaxKi 1274 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV;
884 G4Exception(ss, "em0044", JustWarning, ed); 1275 G4Exception(ss, "em0044", JustWarning, ed);
885 } 1276 }
886 1277
887 //....oooOO0OOooo........oooOO0OOooo........oo 1278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 1279
889 void G4VEmProcess::ProcessDescription(std::ost 1280 void G4VEmProcess::ProcessDescription(std::ostream& out) const
890 { 1281 {
891 if(nullptr != particle) { << 1282 if(particle) {
892 StreamInfo(out, *particle, true); 1283 StreamInfo(out, *particle, true);
893 } 1284 }
894 } 1285 }
895 1286
896 //....oooOO0OOooo........oooOO0OOooo........oo 1287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 1288