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