| 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" 68 #include "G4EmParameters.hh"
69 #include "G4EmProcessSubType.hh" 69 #include "G4EmProcessSubType.hh"
70 #include "G4EmTableUtil.hh" 70 #include "G4EmTableUtil.hh"
71 #include "G4EmUtility.hh" 71 #include "G4EmUtility.hh"
72 #include "G4DNAModelSubType.hh" 72 #include "G4DNAModelSubType.hh"
73 #include "G4GenericIon.hh" 73 #include "G4GenericIon.hh"
74 #include "G4Log.hh" 74 #include "G4Log.hh"
75 #include <iostream> 75 #include <iostream>
76 76
77 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 78
79 G4VEmProcess::G4VEmProcess(const G4String& nam 79 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
80 G4VDiscreteProcess(name, type) 80 G4VDiscreteProcess(name, type)
81 { 81 {
82 theParameters = G4EmParameters::Instance(); 82 theParameters = G4EmParameters::Instance();
83 SetVerboseLevel(1); 83 SetVerboseLevel(1);
84 84
85 // Size of tables 85 // Size of tables
86 minKinEnergy = 0.1*CLHEP::keV; 86 minKinEnergy = 0.1*CLHEP::keV;
87 maxKinEnergy = 100.0*CLHEP::TeV; 87 maxKinEnergy = 100.0*CLHEP::TeV;
88 88
89 // default lambda factor 89 // default lambda factor
90 invLambdaFactor = 1.0/lambdaFactor; 90 invLambdaFactor = 1.0/lambdaFactor;
91 91
92 // particle types 92 // particle types
93 theGamma = G4Gamma::Gamma(); 93 theGamma = G4Gamma::Gamma();
94 theElectron = G4Electron::Electron(); 94 theElectron = G4Electron::Electron();
95 thePositron = G4Positron::Positron(); 95 thePositron = G4Positron::Positron();
96 96
97 pParticleChange = &fParticleChange; 97 pParticleChange = &fParticleChange;
98 fParticleChange.SetSecondaryWeightByProcess( 98 fParticleChange.SetSecondaryWeightByProcess(true);
99 secParticles.reserve(5); 99 secParticles.reserve(5);
100 100
101 modelManager = new G4EmModelManager(); 101 modelManager = new G4EmModelManager();
102 lManager = G4LossTableManager::Instance(); 102 lManager = G4LossTableManager::Instance();
103 lManager->Register(this); 103 lManager->Register(this);
104 isTheMaster = lManager->IsMaster(); 104 isTheMaster = lManager->IsMaster();
105 G4LossTableBuilder* bld = lManager->GetTable 105 G4LossTableBuilder* bld = lManager->GetTableBuilder();
106 theDensityFactor = bld->GetDensityFactors(); 106 theDensityFactor = bld->GetDensityFactors();
107 theDensityIdx = bld->GetCoupleIndexes(); 107 theDensityIdx = bld->GetCoupleIndexes();
108 } 108 }
109 109
110 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 111
112 G4VEmProcess::~G4VEmProcess() 112 G4VEmProcess::~G4VEmProcess()
113 { 113 {
114 if(isTheMaster) { 114 if(isTheMaster) {
115 delete theData; 115 delete theData;
116 delete theEnergyOfCrossSectionMax; 116 delete theEnergyOfCrossSectionMax;
117 } 117 }
118 delete modelManager; 118 delete modelManager;
119 delete biasManager; 119 delete biasManager;
120 lManager->DeRegister(this); 120 lManager->DeRegister(this);
121 } 121 }
122 122
123 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 124
125 void G4VEmProcess::AddEmModel(G4int order, G4V 125 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* ptr,
126 const G4Region* 126 const G4Region* region)
127 { 127 {
128 if(nullptr == ptr) { return; } 128 if(nullptr == ptr) { return; }
129 G4VEmFluctuationModel* fm = nullptr; 129 G4VEmFluctuationModel* fm = nullptr;
130 modelManager->AddEmModel(order, ptr, fm, reg 130 modelManager->AddEmModel(order, ptr, fm, region);
131 ptr->SetParticleChange(pParticleChange); 131 ptr->SetParticleChange(pParticleChange);
132 } 132 }
133 133
134 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 135
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, 136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, G4int)
137 { 137 {
138 if(nullptr == ptr) { return; } 138 if(nullptr == ptr) { return; }
139 if(!emModels.empty()) { 139 if(!emModels.empty()) {
140 for(auto & em : emModels) { if(em == ptr) 140 for(auto & em : emModels) { if(em == ptr) { return; } }
141 } 141 }
142 emModels.push_back(ptr); 142 emModels.push_back(ptr);
143 } 143 }
144 144
145 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 146
147 void G4VEmProcess::PreparePhysicsTable(const G 147 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
148 { 148 {
149 if(nullptr == particle) { SetParticle(&part) 149 if(nullptr == particle) { SetParticle(&part); }
150 150
151 if(part.GetParticleType() == "nucleus" && 151 if(part.GetParticleType() == "nucleus" &&
152 part.GetParticleSubType() == "generic") { 152 part.GetParticleSubType() == "generic") {
153 153
154 G4String pname = part.GetParticleName(); 154 G4String pname = part.GetParticleName();
155 if(pname != "deuteron" && pname != "triton 155 if(pname != "deuteron" && pname != "triton" &&
156 pname != "He3" && pname != "alpha" && p << 156 pname != "alpha" && pname != "alpha+" &&
157 pname != "helium" && pname != "hydrogen 157 pname != "helium" && pname != "hydrogen") {
158 158
159 particle = G4GenericIon::GenericIon(); 159 particle = G4GenericIon::GenericIon();
160 isIon = true; 160 isIon = true;
161 } 161 }
162 } 162 }
163 if(particle != &part) { return; } 163 if(particle != &part) { return; }
164 164
165 lManager->PreparePhysicsTable(&part, this); 165 lManager->PreparePhysicsTable(&part, this);
166 166
167 // for new run 167 // for new run
168 currentCouple = nullptr; 168 currentCouple = nullptr;
169 preStepLambda = 0.0; 169 preStepLambda = 0.0;
170 fLambdaEnergy = 0.0; 170 fLambdaEnergy = 0.0;
171 171
172 InitialiseProcess(particle); 172 InitialiseProcess(particle);
173 173
174 G4LossTableBuilder* bld = lManager->GetTable 174 G4LossTableBuilder* bld = lManager->GetTableBuilder();
175 const G4ProductionCutsTable* theCoupleTable= 175 const G4ProductionCutsTable* theCoupleTable=
176 G4ProductionCutsTable::GetProductionCutsTa 176 G4ProductionCutsTable::GetProductionCutsTable();
177 theCutsGamma = theCoupleTable->GetEnergyC 177 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
178 theCutsElectron = theCoupleTable->GetEnergyC 178 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
179 theCutsPositron = theCoupleTable->GetEnergyC 179 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
180 180
181 // initialisation of the process 181 // initialisation of the process
182 if(!actMinKinEnergy) { minKinEnergy = thePar 182 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
183 if(!actMaxKinEnergy) { maxKinEnergy = thePar 183 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
184 184
185 applyCuts = theParameters->ApplyCuts() 185 applyCuts = theParameters->ApplyCuts();
186 lambdaFactor = theParameters->LambdaFacto 186 lambdaFactor = theParameters->LambdaFactor();
187 invLambdaFactor = 1.0/lambdaFactor; 187 invLambdaFactor = 1.0/lambdaFactor;
188 theParameters->DefineRegParamForEM(this); 188 theParameters->DefineRegParamForEM(this);
189 189
190 // integral option may be disabled 190 // integral option may be disabled
191 if(!theParameters->Integral()) { fXSType = f 191 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; }
192 192
193 // prepare tables 193 // prepare tables
194 if(isTheMaster) { 194 if(isTheMaster) {
195 if(nullptr == theData) { theData = new G4E 195 if(nullptr == theData) { theData = new G4EmDataHandler(2); }
196 196
197 if(buildLambdaTable) { 197 if(buildLambdaTable) {
198 theLambdaTable = theData->MakeTable(0); 198 theLambdaTable = theData->MakeTable(0);
199 bld->InitialiseBaseMaterials(theLambdaTa 199 bld->InitialiseBaseMaterials(theLambdaTable);
200 } 200 }
201 // high energy table 201 // high energy table
202 if(minKinEnergyPrim < maxKinEnergy) { 202 if(minKinEnergyPrim < maxKinEnergy) {
203 theLambdaTablePrim = theData->MakeTable( 203 theLambdaTablePrim = theData->MakeTable(1);
204 bld->InitialiseBaseMaterials(theLambdaTa 204 bld->InitialiseBaseMaterials(theLambdaTablePrim);
205 } 205 }
206 } 206 }
207 // models 207 // models
208 baseMat = bld->GetBaseMaterialFlag(); 208 baseMat = bld->GetBaseMaterialFlag();
209 numberOfModels = modelManager->NumberOfModel 209 numberOfModels = modelManager->NumberOfModels();
210 currentModel = modelManager->GetModel(0); 210 currentModel = modelManager->GetModel(0);
211 if(nullptr != lManager->AtomDeexcitation()) 211 if(nullptr != lManager->AtomDeexcitation()) {
212 modelManager->SetFluoFlag(true); 212 modelManager->SetFluoFlag(true);
213 } 213 }
214 // forced biasing 214 // forced biasing
215 if(nullptr != biasManager) { 215 if(nullptr != biasManager) {
216 biasManager->Initialise(part, GetProcessNa 216 biasManager->Initialise(part, GetProcessName(), verboseLevel);
217 biasFlag = false; 217 biasFlag = false;
218 } 218 }
219 219
220 theCuts = 220 theCuts =
221 G4EmTableUtil::PrepareEmProcess(this, part 221 G4EmTableUtil::PrepareEmProcess(this, particle, secondaryParticle,
222 modelManag 222 modelManager, maxKinEnergy,
223 secID, tri 223 secID, tripletID, mainSecondaries,
224 verboseLev 224 verboseLevel, isTheMaster);
225 } 225 }
226 226
227 //....oooOO0OOooo........oooOO0OOooo........oo 227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 228
229 void G4VEmProcess::BuildPhysicsTable(const G4P 229 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 { 230 {
231 if(nullptr == masterProc) { 231 if(nullptr == masterProc) {
232 if(isTheMaster) { masterProc = this; } 232 if(isTheMaster) { masterProc = this; }
233 else { masterProc = static_cast<const G4VE 233 else { masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess());}
234 } 234 }
235 G4int nModels = modelManager->NumberOfModels 235 G4int nModels = modelManager->NumberOfModels();
236 G4bool isLocked = theParameters->IsPrintLock 236 G4bool isLocked = theParameters->IsPrintLocked();
237 G4bool toBuild = (buildLambdaTable || minKin 237 G4bool toBuild = (buildLambdaTable || minKinEnergyPrim < maxKinEnergy);
238 238
239 G4EmTableUtil::BuildEmProcess(this, masterPr 239 G4EmTableUtil::BuildEmProcess(this, masterProc, particle, &part,
240 nModels, verbo 240 nModels, verboseLevel, isTheMaster,
241 isLocked, toBu 241 isLocked, toBuild, baseMat);
242 } 242 }
243 243
244 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 245
246 void G4VEmProcess::BuildLambdaTable() 246 void G4VEmProcess::BuildLambdaTable()
247 { 247 {
248 G4double scale = theParameters->MaxKinEnergy 248 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
249 G4int nbin = 249 G4int nbin =
250 theParameters->NumberOfBinsPerDecade()*G4l 250 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale));
251 if(actBinning) { nbin = std::max(nbin, nLamb 251 if(actBinning) { nbin = std::max(nbin, nLambdaBins); }
252 scale = nbin/G4Log(scale); 252 scale = nbin/G4Log(scale);
253 253
254 G4LossTableBuilder* bld = lManager->GetTable 254 G4LossTableBuilder* bld = lManager->GetTableBuilder();
255 G4EmTableUtil::BuildLambdaTable(this, partic 255 G4EmTableUtil::BuildLambdaTable(this, particle, modelManager,
256 bld, theLamb 256 bld, theLambdaTable, theLambdaTablePrim,
257 minKinEnergy 257 minKinEnergy, minKinEnergyPrim,
258 maxKinEnergy 258 maxKinEnergy, scale, verboseLevel,
259 startFromNul 259 startFromNull, splineFlag);
260 } 260 }
261 261
262 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 263
264 void G4VEmProcess::StreamInfo(std::ostream& ou 264 void G4VEmProcess::StreamInfo(std::ostream& out,
265 const G4ParticleDefinition& 265 const G4ParticleDefinition& part, G4bool rst) const
266 { 266 {
267 G4String indent = (rst ? " " : ""); 267 G4String indent = (rst ? " " : "");
268 out << std::setprecision(6); 268 out << std::setprecision(6);
269 out << G4endl << indent << GetProcessName() 269 out << G4endl << indent << GetProcessName() << ": ";
270 if (!rst) { 270 if (!rst) {
271 out << " for " << part.GetParticleName(); 271 out << " for " << part.GetParticleName();
272 } 272 }
273 if(fXSType != fEmNoIntegral) { out << " XSt 273 if(fXSType != fEmNoIntegral) { out << " XStype:" << fXSType; }
274 if(applyCuts) { out << " applyCuts:1 "; } 274 if(applyCuts) { out << " applyCuts:1 "; }
275 G4int subtype = GetProcessSubType(); << 275 out << " SubType=" << GetProcessSubType();
276 out << " SubType=" << subtype; << 276 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 << 277 out << " BuildTable=" << buildLambdaTable << G4endl;
284 if(buildLambdaTable) { 278 if(buildLambdaTable) {
285 if(particle == &part) { 279 if(particle == &part) {
286 for(auto & v : *theLambdaTable) { 280 for(auto & v : *theLambdaTable) {
287 if(nullptr != v) { 281 if(nullptr != v) {
288 out << " Lambda table from "; 282 out << " Lambda table from ";
289 G4double emin = v->Energy(0); 283 G4double emin = v->Energy(0);
290 G4double emax = v->GetMaxEnergy(); 284 G4double emax = v->GetMaxEnergy();
291 G4int nbin = G4int(v->GetVectorLengt 285 G4int nbin = G4int(v->GetVectorLength() - 1);
292 if(emin > minKinEnergy) { out << "th 286 if(emin > minKinEnergy) { out << "threshold "; }
293 else { out << G4BestUnit(emin,"Energ 287 else { out << G4BestUnit(emin,"Energy"); }
294 out << " to " 288 out << " to "
295 << G4BestUnit(emax,"Energy") 289 << G4BestUnit(emax,"Energy")
296 << ", " << G4lrint(nbin/std::log 290 << ", " << G4lrint(nbin/std::log10(emax/emin))
297 << " bins/decade, spline: " 291 << " bins/decade, spline: "
298 << splineFlag << G4endl; 292 << splineFlag << G4endl;
299 break; 293 break;
300 } 294 }
301 } 295 }
302 } else { 296 } else {
303 out << " Used Lambda table of " 297 out << " Used Lambda table of "
304 << particle->GetParticleName() << G4endl 298 << particle->GetParticleName() << G4endl;
305 } 299 }
306 } 300 }
307 if(minKinEnergyPrim < maxKinEnergy) { 301 if(minKinEnergyPrim < maxKinEnergy) {
308 if(particle == &part) { 302 if(particle == &part) {
309 for(auto & v : *theLambdaTablePrim) { 303 for(auto & v : *theLambdaTablePrim) {
310 if(nullptr != v) { 304 if(nullptr != v) {
311 out << " LambdaPrime table from 305 out << " LambdaPrime table from "
312 << G4BestUnit(v->Energy(0),"Ener 306 << G4BestUnit(v->Energy(0),"Energy")
313 << " to " 307 << " to "
314 << G4BestUnit(v->GetMaxEnergy(), 308 << G4BestUnit(v->GetMaxEnergy(),"Energy")
315 << " in " << v->GetVectorLength( 309 << " in " << v->GetVectorLength()-1
316 << " bins " << G4endl; 310 << " bins " << G4endl;
317 break; 311 break;
318 } 312 }
319 } 313 }
320 } else { 314 } else {
321 out << " Used LambdaPrime table of 315 out << " Used LambdaPrime table of "
322 << particle->GetParticleName() 316 << particle->GetParticleName() << G4endl;
323 } 317 }
324 } 318 }
325 StreamProcessInfo(out); 319 StreamProcessInfo(out);
326 modelManager->DumpModelList(out, verboseLeve 320 modelManager->DumpModelList(out, verboseLevel);
327 321
328 if(verboseLevel > 2 && buildLambdaTable) { 322 if(verboseLevel > 2 && buildLambdaTable) {
329 out << " LambdaTable address= " << th 323 out << " LambdaTable address= " << theLambdaTable << G4endl;
330 if(theLambdaTable && particle == &part) { 324 if(theLambdaTable && particle == &part) {
331 out << (*theLambdaTable) << G4endl; 325 out << (*theLambdaTable) << G4endl;
332 } 326 }
333 } 327 }
334 } 328 }
335 329
336 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337 331
338 void G4VEmProcess::StartTracking(G4Track* trac 332 void G4VEmProcess::StartTracking(G4Track* track)
339 { 333 {
340 // reset parameters for the new track 334 // reset parameters for the new track
341 currentParticle = track->GetParticleDefiniti 335 currentParticle = track->GetParticleDefinition();
342 theNumberOfInteractionLengthLeft = -1.0; 336 theNumberOfInteractionLengthLeft = -1.0;
343 mfpKinEnergy = DBL_MAX; 337 mfpKinEnergy = DBL_MAX;
344 preStepLambda = 0.0; <<
345 338
346 if(isIon) { massRatio = proton_mass_c2/curre 339 if(isIon) { massRatio = proton_mass_c2/currentParticle->GetPDGMass(); }
347 340
348 // forced biasing only for primary particles 341 // forced biasing only for primary particles
349 if(biasManager) { 342 if(biasManager) {
350 if(0 == track->GetParentID()) { 343 if(0 == track->GetParentID()) {
351 // primary particle 344 // primary particle
352 biasFlag = true; 345 biasFlag = true;
353 biasManager->ResetForcedInteraction(); 346 biasManager->ResetForcedInteraction();
354 } 347 }
355 } 348 }
356 } 349 }
357 350
358 //....oooOO0OOooo........oooOO0OOooo........oo 351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 352
360 G4double G4VEmProcess::PostStepGetPhysicalInte 353 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& tr 354 const G4Track& track,
362 G4double previo 355 G4double previousStepSize,
363 G4ForceCondition* 356 G4ForceCondition* condition)
364 { 357 {
365 *condition = NotForced; 358 *condition = NotForced;
366 G4double x = DBL_MAX; 359 G4double x = DBL_MAX;
367 360
368 DefineMaterial(track.GetMaterialCutsCouple() 361 DefineMaterial(track.GetMaterialCutsCouple());
369 preStepKinEnergy = track.GetKineticEnergy(); << 362 preStepKinEnergy = track.GetKineticEnergy();
>> 363 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
370 const G4double scaledEnergy = preStepKinEner 364 const G4double scaledEnergy = preStepKinEnergy*massRatio;
371 SelectModel(scaledEnergy, currentCoupleIndex 365 SelectModel(scaledEnergy, currentCoupleIndex);
372 /* 366 /*
373 G4cout << "PostStepGetPhysicalInteractionLen 367 G4cout << "PostStepGetPhysicalInteractionLength: idx= " << currentCoupleIndex
374 << " couple: " << currentCouple << G 368 << " couple: " << currentCouple << G4endl;
375 */ 369 */
376 if(!currentModel->IsActive(scaledEnergy)) { 370 if(!currentModel->IsActive(scaledEnergy)) {
377 theNumberOfInteractionLengthLeft = -1.0; 371 theNumberOfInteractionLengthLeft = -1.0;
378 currentInteractionLength = DBL_MAX; 372 currentInteractionLength = DBL_MAX;
379 mfpKinEnergy = DBL_MAX; <<
380 preStepLambda = 0.0; <<
381 return x; 373 return x;
382 } 374 }
383 375
384 // forced biasing only for primary particles 376 // forced biasing only for primary particles
385 if(biasManager) { 377 if(biasManager) {
386 if(0 == track.GetParentID()) { 378 if(0 == track.GetParentID()) {
387 if(biasFlag && 379 if(biasFlag &&
388 biasManager->ForcedInteractionRegion( 380 biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
389 return biasManager->GetStepLimit((G4in 381 return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize);
390 } 382 }
391 } 383 }
392 } 384 }
393 385
394 // compute mean free path 386 // compute mean free path
395 << 387 ComputeIntegralLambda(preStepKinEnergy, preStepLogKinEnergy);
396 ComputeIntegralLambda(preStepKinEnergy, trac <<
397 388
398 // zero cross section 389 // zero cross section
399 if(preStepLambda <= 0.0) { 390 if(preStepLambda <= 0.0) {
400 theNumberOfInteractionLengthLeft = -1.0; 391 theNumberOfInteractionLengthLeft = -1.0;
401 currentInteractionLength = DBL_MAX; 392 currentInteractionLength = DBL_MAX;
402 393
403 } else { 394 } else {
404 395
405 // non-zero cross section 396 // non-zero cross section
406 if (theNumberOfInteractionLengthLeft < 0.0 397 if (theNumberOfInteractionLengthLeft < 0.0) {
407 398
408 // beggining of tracking (or just after 399 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Lo 400 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
410 theInitialNumberOfInteractionLength = th 401 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
411 402
412 } else { 403 } else {
413 404
414 theNumberOfInteractionLengthLeft -= 405 theNumberOfInteractionLengthLeft -=
415 previousStepSize/currentInteractionLen 406 previousStepSize/currentInteractionLength;
416 theNumberOfInteractionLengthLeft = 407 theNumberOfInteractionLengthLeft =
417 std::max(theNumberOfInteractionLengthL 408 std::max(theNumberOfInteractionLengthLeft, 0.0);
418 } 409 }
419 410
420 // new mean free path and step limit for t 411 // new mean free path and step limit for the next step
421 currentInteractionLength = 1.0/preStepLamb 412 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * cur 413 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 } 414 }
424 return x; 415 return x;
425 } 416 }
426 417
427 //....oooOO0OOooo........oooOO0OOooo........oo 418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428 419
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 420 void G4VEmProcess::ComputeIntegralLambda(G4double e, G4double loge)
430 { 421 {
431 if (fXSType == fEmNoIntegral) { << 422 if(fXSType == fEmNoIntegral) {
432 preStepLambda = GetCurrentLambda(e, LogEki << 423 preStepLambda = GetCurrentLambda(e, loge);
433 424
434 } else if (fXSType == fEmIncreasing) { << 425 } else if(fXSType == fEmIncreasing) {
435 if(e*invLambdaFactor < mfpKinEnergy) { 426 if(e*invLambdaFactor < mfpKinEnergy) {
436 preStepLambda = GetCurrentLambda(e, LogE << 427 mfpKinEnergy = e;
437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 428 preStepLambda = GetCurrentLambda(e, loge);
438 } 429 }
439 430
440 } else if(fXSType == fEmDecreasing) { 431 } else if(fXSType == fEmDecreasing) {
441 if(e < mfpKinEnergy) { 432 if(e < mfpKinEnergy) {
442 const G4double e1 = e*lambdaFactor; 433 const G4double e1 = e*lambdaFactor;
443 preStepLambda = GetCurrentLambda(e1); 434 preStepLambda = GetCurrentLambda(e1);
444 mfpKinEnergy = e1; 435 mfpKinEnergy = e1;
445 } 436 }
446 437
447 } else if(fXSType == fEmOnePeak) { 438 } else if(fXSType == fEmOnePeak) {
448 const G4double epeak = (*theEnergyOfCrossS 439 const G4double epeak = (*theEnergyOfCrossSectionMax)[currentCoupleIndex];
449 if(e <= epeak) { 440 if(e <= epeak) {
450 if(e*invLambdaFactor < mfpKinEnergy) { 441 if(e*invLambdaFactor < mfpKinEnergy) {
451 preStepLambda = GetCurrentLambda(e, Lo << 442 mfpKinEnergy = e;
452 mfpKinEnergy = (preStepLambda > 0.0) ? << 443 preStepLambda = GetCurrentLambda(e, loge);
453 } 444 }
454 } else if(e < mfpKinEnergy) { 445 } else if(e < mfpKinEnergy) {
455 const G4double e1 = std::max(epeak, e*la 446 const G4double e1 = std::max(epeak, e*lambdaFactor);
456 preStepLambda = GetCurrentLambda(e1); << 447 preStepLambda = GetCurrentLambda(e1);
457 mfpKinEnergy = e1; 448 mfpKinEnergy = e1;
458 } 449 }
>> 450
459 } else { 451 } else {
460 preStepLambda = GetCurrentLambda(e, LogEki << 452 preStepLambda = GetCurrentLambda(e, loge);
461 } 453 }
462 } 454 }
463 455
464 //....oooOO0OOooo........oooOO0OOooo........oo 456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 457
466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 458 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 459 const G4Step& step)
468 { 460 {
469 // clear number of interaction lengths in an 461 // clear number of interaction lengths in any case
470 theNumberOfInteractionLengthLeft = -1.0; 462 theNumberOfInteractionLengthLeft = -1.0;
471 mfpKinEnergy = DBL_MAX; 463 mfpKinEnergy = DBL_MAX;
472 464
473 fParticleChange.InitializeForPostStep(track) 465 fParticleChange.InitializeForPostStep(track);
474 466
475 // Do not make anything if particle is stopp 467 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt! 468 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) 469 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478 470
479 const G4double finalT = track.GetKineticEner 471 const G4double finalT = track.GetKineticEnergy();
480 472
481 // forced process - should happen only once 473 // forced process - should happen only once per track
482 if(biasFlag) { 474 if(biasFlag) {
483 if(biasManager->ForcedInteractionRegion((G 475 if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
484 biasFlag = false; 476 biasFlag = false;
485 } 477 }
486 } 478 }
487 479
488 // check active and select model 480 // check active and select model
489 const G4double scaledEnergy = finalT*massRat 481 const G4double scaledEnergy = finalT*massRatio;
490 SelectModel(scaledEnergy, currentCoupleIndex 482 SelectModel(scaledEnergy, currentCoupleIndex);
491 if(!currentModel->IsActive(scaledEnergy)) { 483 if(!currentModel->IsActive(scaledEnergy)) { return &fParticleChange; }
492 484
493 // Integral approach 485 // Integral approach
494 if (fXSType != fEmNoIntegral) { 486 if (fXSType != fEmNoIntegral) {
495 const G4double logFinalT = 487 const G4double logFinalT =
496 track.GetDynamicParticle()->GetLogKineti 488 track.GetDynamicParticle()->GetLogKineticEnergy();
497 const G4double lx = std::max(GetCurrentLam 489 const G4double lx = std::max(GetCurrentLambda(finalT, logFinalT), 0.0);
498 #ifdef G4VERBOSE 490 #ifdef G4VERBOSE
499 if(preStepLambda < lx && 1 < verboseLevel) 491 if(preStepLambda < lx && 1 < verboseLevel) {
500 G4cout << "WARNING: for " << currentPart 492 G4cout << "WARNING: for " << currentParticle->GetParticleName()
501 << " and " << GetProcessName() << 493 << " and " << GetProcessName() << " E(MeV)= " << finalT/MeV
502 << " preLambda= " << preStepLambd 494 << " preLambda= " << preStepLambda
503 << " < " << lx << " (postLambda) 495 << " < " << lx << " (postLambda) " << G4endl;
504 } 496 }
505 #endif 497 #endif
506 // if false interaction then use new cross 498 // if false interaction then use new cross section value
507 // if both values are zero - no interactio 499 // if both values are zero - no interaction
508 if(preStepLambda*G4UniformRand() >= lx) { 500 if(preStepLambda*G4UniformRand() >= lx) {
509 return &fParticleChange; 501 return &fParticleChange;
510 } 502 }
511 } 503 }
512 504
513 // define new weight for primary and seconda 505 // define new weight for primary and secondaries
514 G4double weight = fParticleChange.GetParentW 506 G4double weight = fParticleChange.GetParentWeight();
515 if(weightFlag) { 507 if(weightFlag) {
516 weight /= biasFactor; 508 weight /= biasFactor;
517 fParticleChange.ProposeWeight(weight); 509 fParticleChange.ProposeWeight(weight);
518 } 510 }
519 511
520 #ifdef G4VERBOSE 512 #ifdef G4VERBOSE
521 if(1 < verboseLevel) { 513 if(1 < verboseLevel) {
522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 514 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV 515 << finalT/MeV
524 << " MeV; model= (" << currentModel 516 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnerg 517 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl; 518 << G4endl;
527 } 519 }
528 #endif 520 #endif
529 521
530 // sample secondaries 522 // sample secondaries
531 secParticles.clear(); 523 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticle 524 currentModel->SampleSecondaries(&secParticles,
533 currentCoupl 525 currentCouple,
534 track.GetDyn 526 track.GetDynamicParticle(),
535 (*theCuts)[c 527 (*theCuts)[currentCoupleIndex]);
536 528
537 G4int num0 = (G4int)secParticles.size(); 529 G4int num0 = (G4int)secParticles.size();
538 530
539 // splitting or Russian roulette 531 // splitting or Russian roulette
540 if(biasManager) { 532 if(biasManager) {
541 if(biasManager->SecondaryBiasingRegion((G4 533 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
542 G4double eloss = 0.0; 534 G4double eloss = 0.0;
543 weight *= biasManager->ApplySecondaryBia 535 weight *= biasManager->ApplySecondaryBiasing(
544 secParticles, track, currentModel, &fP 536 secParticles, track, currentModel, &fParticleChange, eloss,
545 (G4int)currentCoupleIndex, (*theCuts)[ 537 (G4int)currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546 step.GetPostStepPoint()->GetSafety()); 538 step.GetPostStepPoint()->GetSafety());
547 if(eloss > 0.0) { 539 if(eloss > 0.0) {
548 eloss += fParticleChange.GetLocalEnerg 540 eloss += fParticleChange.GetLocalEnergyDeposit();
549 fParticleChange.ProposeLocalEnergyDepo 541 fParticleChange.ProposeLocalEnergyDeposit(eloss);
550 } 542 }
551 } 543 }
552 } 544 }
553 545
554 // save secondaries 546 // save secondaries
555 G4int num = (G4int)secParticles.size(); 547 G4int num = (G4int)secParticles.size();
556 if(num > 0) { 548 if(num > 0) {
557 549
558 fParticleChange.SetNumberOfSecondaries(num 550 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEn 551 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime(); 552 G4double time = track.GetGlobalTime();
561 553
562 G4int n1(0), n2(0); 554 G4int n1(0), n2(0);
563 if(num0 > mainSecondaries) { << 555 if(num > mainSecondaries) {
564 currentModel->FillNumberOfSecondaries(n1 556 currentModel->FillNumberOfSecondaries(n1, n2);
565 } 557 }
566 558
567 for (G4int i=0; i<num; ++i) { 559 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i]; 560 G4DynamicParticle* dp = secParticles[i];
569 if (nullptr != dp) { 561 if (nullptr != dp) {
570 const G4ParticleDefinition* p = dp->Ge 562 const G4ParticleDefinition* p = dp->GetParticleDefinition();
571 G4double e = dp->GetKineticEnergy(); 563 G4double e = dp->GetKineticEnergy();
572 G4bool good = true; 564 G4bool good = true;
573 if(applyCuts) { 565 if(applyCuts) {
574 if (p == theGamma) { 566 if (p == theGamma) {
575 if (e < (*theCutsGamma)[currentCou 567 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576 568
577 } else if (p == theElectron) { 569 } else if (p == theElectron) {
578 if (e < (*theCutsElectron)[current 570 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579 571
580 } else if (p == thePositron) { 572 } else if (p == thePositron) {
581 if (electron_mass_c2 < (*theCutsGa 573 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582 e < (*theCutsPositron)[current 574 e < (*theCutsPositron)[currentCoupleIndex]) {
583 good = false; 575 good = false;
584 e += 2.0*electron_mass_c2; 576 e += 2.0*electron_mass_c2;
585 } 577 }
586 } 578 }
587 // added secondary if it is good 579 // added secondary if it is good
588 } 580 }
589 if (good) { 581 if (good) {
590 G4Track* t = new G4Track(dp, time, t 582 G4Track* t = new G4Track(dp, time, track.GetPosition());
591 t->SetTouchableHandle(track.GetTouch 583 t->SetTouchableHandle(track.GetTouchableHandle());
592 if (biasManager) { 584 if (biasManager) {
593 t->SetWeight(weight * biasManager- 585 t->SetWeight(weight * biasManager->GetWeight(i));
594 } else { 586 } else {
595 t->SetWeight(weight); 587 t->SetWeight(weight);
596 } 588 }
597 pParticleChange->AddSecondary(t); 589 pParticleChange->AddSecondary(t);
598 590
599 // define type of secondary 591 // define type of secondary
600 if(i < mainSecondaries) { 592 if(i < mainSecondaries) {
601 t->SetCreatorModelID(secID); 593 t->SetCreatorModelID(secID);
602 if(GetProcessSubType() == fCompton 594 if(GetProcessSubType() == fComptonScattering && p == theGamma) {
603 t->SetCreatorModelID(_ComptonGam 595 t->SetCreatorModelID(_ComptonGamma);
604 } 596 }
605 } else if(i < mainSecondaries + n1) 597 } else if(i < mainSecondaries + n1) {
606 t->SetCreatorModelID(tripletID); 598 t->SetCreatorModelID(tripletID);
607 } else if(i < mainSecondaries + n1 + 599 } else if(i < mainSecondaries + n1 + n2) {
608 t->SetCreatorModelID(_IonRecoil); 600 t->SetCreatorModelID(_IonRecoil);
609 } else { 601 } else {
610 if(i < num0) { 602 if(i < num0) {
611 if(p == theGamma) { 603 if(p == theGamma) {
612 t->SetCreatorModelID(fluoID); 604 t->SetCreatorModelID(fluoID);
613 } else { 605 } else {
614 t->SetCreatorModelID(augerID); 606 t->SetCreatorModelID(augerID);
615 } 607 }
616 } else { 608 } else {
617 t->SetCreatorModelID(biasID); << 609 t->SetCreatorModelID(secID);
618 } 610 }
619 } 611 }
620 /* 612 /*
621 G4cout << "Secondary(post step) has 613 G4cout << "Secondary(post step) has weight " << t->GetWeight()
622 << ", Ekin= " << t->GetKineti 614 << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV "
623 << GetProcessName() << " fluo 615 << GetProcessName() << " fluoID= " << fluoID
624 << " augerID= " << augerID << 616 << " augerID= " << augerID <<G4endl;
625 */ 617 */
626 } else { 618 } else {
627 delete dp; 619 delete dp;
628 edep += e; 620 edep += e;
629 } 621 }
630 } 622 }
631 } 623 }
632 fParticleChange.ProposeLocalEnergyDeposit( 624 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 } 625 }
634 626
635 if(0.0 == fParticleChange.GetProposedKinetic 627 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636 fAlive == fParticleChange.GetTrackStatus( 628 fAlive == fParticleChange.GetTrackStatus()) {
637 if(particle->GetProcessManager()->GetAtRes 629 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638 { fParticleChange.ProposeTrackStatus( 630 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639 else { fParticleChange.ProposeTrackStatus( 631 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640 } 632 }
641 633
642 return &fParticleChange; 634 return &fParticleChange;
643 } 635 }
644 636
645 //....oooOO0OOooo........oooOO0OOooo........oo 637 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 638
647 G4bool G4VEmProcess::StorePhysicsTable(const G 639 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G 640 const G4String& directory,
649 G4bool 641 G4bool ascii)
650 { 642 {
651 if(!isTheMaster || part != particle) { retur 643 if(!isTheMaster || part != particle) { return true; }
652 if(G4EmTableUtil::StoreTable(this, part, the 644 if(G4EmTableUtil::StoreTable(this, part, theLambdaTable,
653 directory, "Lambda", 645 directory, "Lambda",
654 verboseLevel, a 646 verboseLevel, ascii) &&
655 G4EmTableUtil::StoreTable(this, part, the 647 G4EmTableUtil::StoreTable(this, part, theLambdaTablePrim,
656 directory, "LambdaPrim", 648 directory, "LambdaPrim",
657 verboseLevel, a 649 verboseLevel, ascii)) {
658 return true; 650 return true;
659 } 651 }
660 return false; 652 return false;
661 } 653 }
662 654
663 //....oooOO0OOooo........oooOO0OOooo........oo 655 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664 656
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 657 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 cons 658 const G4String& dir,
667 G4bo 659 G4bool ascii)
668 { 660 {
669 if(!isTheMaster || part != particle) { retur 661 if(!isTheMaster || part != particle) { return true; }
670 G4bool yes = true; 662 G4bool yes = true;
671 if(buildLambdaTable) { 663 if(buildLambdaTable) {
672 yes = G4EmTableUtil::RetrieveTable(this, p 664 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTable, dir,
673 "Lambda 665 "Lambda", verboseLevel,
674 ascii, 666 ascii, splineFlag);
675 } 667 }
676 if(yes && minKinEnergyPrim < maxKinEnergy) { 668 if(yes && minKinEnergyPrim < maxKinEnergy) {
677 yes = G4EmTableUtil::RetrieveTable(this, p 669 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTablePrim, dir,
678 "Lambda 670 "LambdaPrim", verboseLevel,
679 ascii, 671 ascii, splineFlag);
680 } 672 }
681 return yes; 673 return yes;
682 } 674 }
683 675
684 //....oooOO0OOooo........oooOO0OOooo........oo 676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685 677
686 G4double G4VEmProcess::GetCrossSection(G4doubl 678 G4double G4VEmProcess::GetCrossSection(G4double kinEnergy,
687 const G 679 const G4MaterialCutsCouple* couple)
688 { 680 {
689 CurrentSetup(couple, kinEnergy); 681 CurrentSetup(couple, kinEnergy);
690 return GetCurrentLambda(kinEnergy, G4Log(kin 682 return GetCurrentLambda(kinEnergy, G4Log(kinEnergy));
691 } 683 }
692 684
693 //....oooOO0OOooo........oooOO0OOooo........oo 685 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694 686
695 G4double G4VEmProcess::GetMeanFreePath(const G 687 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696 G4doubl 688 G4double,
697 G4Force 689 G4ForceCondition* condition)
698 { 690 {
699 *condition = NotForced; 691 *condition = NotForced;
700 return G4VEmProcess::MeanFreePath(track); 692 return G4VEmProcess::MeanFreePath(track);
701 } 693 }
702 694
703 //....oooOO0OOooo........oooOO0OOooo........oo 695 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 696
705 G4double 697 G4double
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou 698 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kinEnergy,
707 G4dou 699 G4double Z, G4double A, G4double cut)
708 { 700 {
709 SelectModel(kinEnergy, currentCoupleIndex); 701 SelectModel(kinEnergy, currentCoupleIndex);
710 return (currentModel) ? 702 return (currentModel) ?
711 currentModel->ComputeCrossSectionPerAtom(c 703 currentModel->ComputeCrossSectionPerAtom(currentParticle, kinEnergy,
712 Z 704 Z, A, cut) : 0.0;
713 } 705 }
714 706
715 //....oooOO0OOooo........oooOO0OOooo........oo 707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 708
717 G4PhysicsVector* 709 G4PhysicsVector*
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate 710 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple)
719 { 711 {
720 DefineMaterial(couple); 712 DefineMaterial(couple);
721 G4PhysicsVector* newv = new G4PhysicsLogVect 713 G4PhysicsVector* newv = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy,
722 714 nLambdaBins, splineFlag);
723 return newv; 715 return newv;
724 } 716 }
725 717
726 //....oooOO0OOooo........oooOO0OOooo........oo 718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
727 719
728 const G4Element* G4VEmProcess::GetCurrentEleme 720 const G4Element* G4VEmProcess::GetCurrentElement() const
729 { 721 {
730 return (nullptr != currentModel) ? 722 return (nullptr != currentModel) ?
731 currentModel->GetCurrentElement(currentMat 723 currentModel->GetCurrentElement(currentMaterial) : nullptr;
732 } 724 }
733 725
734 //....oooOO0OOooo........oooOO0OOooo........oo 726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735 727
736 const G4Element* G4VEmProcess::GetTargetElemen 728 const G4Element* G4VEmProcess::GetTargetElement() const
737 { 729 {
738 return (nullptr != currentModel) ? 730 return (nullptr != currentModel) ?
739 currentModel->GetCurrentElement(currentMat 731 currentModel->GetCurrentElement(currentMaterial) : nullptr;
740 } 732 }
741 733
742 //....oooOO0OOooo........oooOO0OOooo........oo 734 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743 735
744 const G4Isotope* G4VEmProcess::GetTargetIsotop 736 const G4Isotope* G4VEmProcess::GetTargetIsotope() const
745 { 737 {
746 return (nullptr != currentModel) ? 738 return (nullptr != currentModel) ?
747 currentModel->GetCurrentIsotope(GetCurrent 739 currentModel->GetCurrentIsotope(GetCurrentElement()) : nullptr;
748 } 740 }
749 741
750 //....oooOO0OOooo........oooOO0OOooo........oo 742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 743
752 void G4VEmProcess::SetCrossSectionBiasingFacto 744 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 { 745 {
754 if(f > 0.0) { 746 if(f > 0.0) {
755 biasFactor = f; 747 biasFactor = f;
756 weightFlag = flag; 748 weightFlag = flag;
757 if(1 < verboseLevel) { 749 if(1 < verboseLevel) {
758 G4cout << "### SetCrossSectionBiasingFac 750 G4cout << "### SetCrossSectionBiasingFactor: for "
759 << particle->GetParticleName() 751 << particle->GetParticleName()
760 << " and process " << GetProcessN 752 << " and process " << GetProcessName()
761 << " biasFactor= " << f << " weig 753 << " biasFactor= " << f << " weightFlag= " << flag
762 << G4endl; 754 << G4endl;
763 } 755 }
764 } 756 }
765 } 757 }
766 758
767 //....oooOO0OOooo........oooOO0OOooo........oo 759 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768 760
769 void 761 void
770 G4VEmProcess::ActivateForcedInteraction(G4doub 762 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771 G4bool 763 G4bool flag)
772 { 764 {
773 if(nullptr == biasManager) { biasManager = n 765 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
774 if(1 < verboseLevel) { 766 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: 767 G4cout << "### ActivateForcedInteraction: for "
776 << particle->GetParticleName() 768 << particle->GetParticleName()
777 << " and process " << GetProcessNam 769 << " and process " << GetProcessName()
778 << " length(mm)= " << length/mm 770 << " length(mm)= " << length/mm
779 << " in G4Region <" << r 771 << " in G4Region <" << r
780 << "> weightFlag= " << flag 772 << "> weightFlag= " << flag
781 << G4endl; 773 << G4endl;
782 } 774 }
783 weightFlag = flag; 775 weightFlag = flag;
784 biasManager->ActivateForcedInteraction(lengt 776 biasManager->ActivateForcedInteraction(length, r);
785 } 777 }
786 778
787 //....oooOO0OOooo........oooOO0OOooo........oo 779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788 780
789 void 781 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G 782 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791 G4double factor, 783 G4double factor,
792 G4double energyLimit) 784 G4double energyLimit)
793 { 785 {
794 if (0.0 <= factor) { 786 if (0.0 <= factor) {
795 787
796 // Range cut can be applied only for e- 788 // Range cut can be applied only for e-
797 if(0.0 == factor && secondaryParticle != G 789 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798 { return; } 790 { return; }
799 791
800 if(!biasManager) { biasManager = new G4EmB 792 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801 biasManager->ActivateSecondaryBiasing(regi 793 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802 if(1 < verboseLevel) { 794 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: 795 G4cout << "### ActivateSecondaryBiasing: for "
804 << " process " << GetProcessName() 796 << " process " << GetProcessName()
805 << " factor= " << factor 797 << " factor= " << factor
806 << " in G4Region <" << region 798 << " in G4Region <" << region
807 << "> energyLimit(MeV)= " << energyLimi 799 << "> energyLimit(MeV)= " << energyLimit/MeV
808 << G4endl; 800 << G4endl;
809 } 801 }
810 } 802 }
811 } 803 }
812 804
813 //....oooOO0OOooo........oooOO0OOooo........oo 805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 806
815 void G4VEmProcess::SetLambdaBinning(G4int n) 807 void G4VEmProcess::SetLambdaBinning(G4int n)
816 { 808 {
817 if(5 < n && n < 10000000) { 809 if(5 < n && n < 10000000) {
818 nLambdaBins = n; 810 nLambdaBins = n;
819 actBinning = true; 811 actBinning = true;
820 } else { 812 } else {
821 G4double e = (G4double)n; 813 G4double e = (G4double)n;
822 PrintWarning("SetLambdaBinning", e); 814 PrintWarning("SetLambdaBinning", e);
823 } 815 }
824 } 816 }
825 817
826 //....oooOO0OOooo........oooOO0OOooo........oo 818 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827 819
828 void G4VEmProcess::SetMinKinEnergy(G4double e) 820 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 { 821 {
830 if(1.e-3*eV < e && e < maxKinEnergy) { 822 if(1.e-3*eV < e && e < maxKinEnergy) {
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma 823 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e)
832 /G4Log(maxKinEnergy/ 824 /G4Log(maxKinEnergy/minKinEnergy));
833 minKinEnergy = e; 825 minKinEnergy = e;
834 actMinKinEnergy = true; 826 actMinKinEnergy = true;
835 } else { PrintWarning("SetMinKinEnergy", e); 827 } else { PrintWarning("SetMinKinEnergy", e); }
836 } 828 }
837 829
838 //....oooOO0OOooo........oooOO0OOooo........oo 830 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839 831
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) 832 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 { 833 {
842 if(minKinEnergy < e && e < 1.e+6*TeV) { 834 if(minKinEnergy < e && e < 1.e+6*TeV) {
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ 835 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy)
844 /G4Log(maxKinEnergy/ 836 /G4Log(maxKinEnergy/minKinEnergy));
845 maxKinEnergy = e; 837 maxKinEnergy = e;
846 actMaxKinEnergy = true; 838 actMaxKinEnergy = true;
847 } else { PrintWarning("SetMaxKinEnergy", e); 839 } else { PrintWarning("SetMaxKinEnergy", e); }
848 } 840 }
849 841
850 //....oooOO0OOooo........oooOO0OOooo........oo 842 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851 843
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl 844 void G4VEmProcess::SetMinKinEnergyPrim(G4double e)
853 { 845 {
854 if(theParameters->MinKinEnergy() <= e && 846 if(theParameters->MinKinEnergy() <= e &&
855 e <= theParameters->MaxKinEnergy()) { min 847 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; }
856 else { PrintWarning("SetMinKinEnergyPrim", e 848 else { PrintWarning("SetMinKinEnergyPrim", e); }
857 } 849 }
858 850
859 //....oooOO0OOooo........oooOO0OOooo........oo 851 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
860 852
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const 853 G4VEmProcess* G4VEmProcess::GetEmProcess(const G4String& nam)
862 { 854 {
863 return (nam == GetProcessName()) ? this : nu 855 return (nam == GetProcessName()) ? this : nullptr;
864 } 856 }
865 857
866 //....oooOO0OOooo........oooOO0OOooo........oo 858 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
867 859
868 G4double G4VEmProcess::PolarAngleLimit() const 860 G4double G4VEmProcess::PolarAngleLimit() const
869 { 861 {
870 return theParameters->MscThetaLimit(); 862 return theParameters->MscThetaLimit();
871 } 863 }
872 864
873 //....oooOO0OOooo........oooOO0OOooo........oo 865 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
874 866
875 void G4VEmProcess::PrintWarning(G4String tit, 867 void G4VEmProcess::PrintWarning(G4String tit, G4double val)
876 { 868 {
877 G4String ss = "G4VEmProcess::" + tit; 869 G4String ss = "G4VEmProcess::" + tit;
878 G4ExceptionDescription ed; 870 G4ExceptionDescription ed;
879 ed << "Parameter is out of range: " << val 871 ed << "Parameter is out of range: " << val
880 << " it will have no effect!\n" << " Pro 872 << " it will have no effect!\n" << " Process "
881 << GetProcessName() << " nbins= " << the 873 << GetProcessName() << " nbins= " << theParameters->NumberOfBins()
882 << " Emin(keV)= " << theParameters->MinKi 874 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV
883 << " Emax(GeV)= " << theParameters->MaxKi 875 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV;
884 G4Exception(ss, "em0044", JustWarning, ed); 876 G4Exception(ss, "em0044", JustWarning, ed);
885 } 877 }
886 878
887 //....oooOO0OOooo........oooOO0OOooo........oo 879 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 880
889 void G4VEmProcess::ProcessDescription(std::ost 881 void G4VEmProcess::ProcessDescription(std::ostream& out) const
890 { 882 {
891 if(nullptr != particle) { 883 if(nullptr != particle) {
892 StreamInfo(out, *particle, true); 884 StreamInfo(out, *particle, true);
893 } 885 }
894 } 886 }
895 887
896 //....oooOO0OOooo........oooOO0OOooo........oo 888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 889