| 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; 338 preStepLambda = 0.0;
345 339
346 if(isIon) { massRatio = proton_mass_c2/curre 340 if(isIon) { massRatio = proton_mass_c2/currentParticle->GetPDGMass(); }
347 341
348 // forced biasing only for primary particles 342 // forced biasing only for primary particles
349 if(biasManager) { 343 if(biasManager) {
350 if(0 == track->GetParentID()) { 344 if(0 == track->GetParentID()) {
351 // primary particle 345 // primary particle
352 biasFlag = true; 346 biasFlag = true;
353 biasManager->ResetForcedInteraction(); 347 biasManager->ResetForcedInteraction();
354 } 348 }
355 } 349 }
356 } 350 }
357 351
358 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 353
360 G4double G4VEmProcess::PostStepGetPhysicalInte 354 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& tr 355 const G4Track& track,
362 G4double previo 356 G4double previousStepSize,
363 G4ForceCondition* 357 G4ForceCondition* condition)
364 { 358 {
365 *condition = NotForced; 359 *condition = NotForced;
366 G4double x = DBL_MAX; 360 G4double x = DBL_MAX;
367 361
368 DefineMaterial(track.GetMaterialCutsCouple() 362 DefineMaterial(track.GetMaterialCutsCouple());
369 preStepKinEnergy = track.GetKineticEnergy(); 363 preStepKinEnergy = track.GetKineticEnergy();
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; 373 mfpKinEnergy = DBL_MAX;
380 preStepLambda = 0.0; 374 preStepLambda = 0.0;
381 return x; 375 return x;
382 } 376 }
383 377
384 // forced biasing only for primary particles 378 // forced biasing only for primary particles
385 if(biasManager) { 379 if(biasManager) {
386 if(0 == track.GetParentID()) { 380 if(0 == track.GetParentID()) {
387 if(biasFlag && 381 if(biasFlag &&
388 biasManager->ForcedInteractionRegion( 382 biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
389 return biasManager->GetStepLimit((G4in 383 return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize);
390 } 384 }
391 } 385 }
392 } 386 }
393 387
394 // compute mean free path 388 // compute mean free path
395 389
396 ComputeIntegralLambda(preStepKinEnergy, trac 390 ComputeIntegralLambda(preStepKinEnergy, track);
397 391
398 // zero cross section 392 // zero cross section
399 if(preStepLambda <= 0.0) { 393 if(preStepLambda <= 0.0) {
400 theNumberOfInteractionLengthLeft = -1.0; 394 theNumberOfInteractionLengthLeft = -1.0;
401 currentInteractionLength = DBL_MAX; 395 currentInteractionLength = DBL_MAX;
402 396
403 } else { 397 } else {
404 398
405 // non-zero cross section 399 // non-zero cross section
406 if (theNumberOfInteractionLengthLeft < 0.0 400 if (theNumberOfInteractionLengthLeft < 0.0) {
407 401
408 // beggining of tracking (or just after 402 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Lo 403 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
410 theInitialNumberOfInteractionLength = th 404 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
411 405
412 } else { 406 } else {
413 407
414 theNumberOfInteractionLengthLeft -= 408 theNumberOfInteractionLengthLeft -=
415 previousStepSize/currentInteractionLen 409 previousStepSize/currentInteractionLength;
416 theNumberOfInteractionLengthLeft = 410 theNumberOfInteractionLengthLeft =
417 std::max(theNumberOfInteractionLengthL 411 std::max(theNumberOfInteractionLengthLeft, 0.0);
418 } 412 }
419 413
420 // new mean free path and step limit for t 414 // new mean free path and step limit for the next step
421 currentInteractionLength = 1.0/preStepLamb 415 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * cur 416 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 } 417 }
424 return x; 418 return x;
425 } 419 }
426 420
427 //....oooOO0OOooo........oooOO0OOooo........oo 421 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428 422
429 void G4VEmProcess::ComputeIntegralLambda(G4dou 423 void G4VEmProcess::ComputeIntegralLambda(G4double e, const G4Track& track)
430 { 424 {
431 if (fXSType == fEmNoIntegral) { 425 if (fXSType == fEmNoIntegral) {
432 preStepLambda = GetCurrentLambda(e, LogEki 426 preStepLambda = GetCurrentLambda(e, LogEkin(track));
433 427
434 } else if (fXSType == fEmIncreasing) { 428 } else if (fXSType == fEmIncreasing) {
435 if(e*invLambdaFactor < mfpKinEnergy) { 429 if(e*invLambdaFactor < mfpKinEnergy) {
436 preStepLambda = GetCurrentLambda(e, LogE 430 preStepLambda = GetCurrentLambda(e, LogEkin(track));
437 mfpKinEnergy = (preStepLambda > 0.0) ? e 431 mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
438 } 432 }
439 433
440 } else if(fXSType == fEmDecreasing) { 434 } else if(fXSType == fEmDecreasing) {
441 if(e < mfpKinEnergy) { 435 if(e < mfpKinEnergy) {
442 const G4double e1 = e*lambdaFactor; 436 const G4double e1 = e*lambdaFactor;
443 preStepLambda = GetCurrentLambda(e1); 437 preStepLambda = GetCurrentLambda(e1);
444 mfpKinEnergy = e1; 438 mfpKinEnergy = e1;
445 } 439 }
446 440
447 } else if(fXSType == fEmOnePeak) { 441 } else if(fXSType == fEmOnePeak) {
448 const G4double epeak = (*theEnergyOfCrossS 442 const G4double epeak = (*theEnergyOfCrossSectionMax)[currentCoupleIndex];
449 if(e <= epeak) { 443 if(e <= epeak) {
450 if(e*invLambdaFactor < mfpKinEnergy) { 444 if(e*invLambdaFactor < mfpKinEnergy) {
451 preStepLambda = GetCurrentLambda(e, Lo 445 preStepLambda = GetCurrentLambda(e, LogEkin(track));
452 mfpKinEnergy = (preStepLambda > 0.0) ? 446 mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
453 } 447 }
454 } else if(e < mfpKinEnergy) { 448 } else if(e < mfpKinEnergy) {
455 const G4double e1 = std::max(epeak, e*la 449 const G4double e1 = std::max(epeak, e*lambdaFactor);
456 preStepLambda = GetCurrentLambda(e1); 450 preStepLambda = GetCurrentLambda(e1);
457 mfpKinEnergy = e1; 451 mfpKinEnergy = e1;
458 } 452 }
459 } else { 453 } else {
460 preStepLambda = GetCurrentLambda(e, LogEki 454 preStepLambda = GetCurrentLambda(e, LogEkin(track));
461 } 455 }
462 } 456 }
463 457
464 //....oooOO0OOooo........oooOO0OOooo........oo 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 459
466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 460 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 461 const G4Step& step)
468 { 462 {
469 // clear number of interaction lengths in an 463 // clear number of interaction lengths in any case
470 theNumberOfInteractionLengthLeft = -1.0; 464 theNumberOfInteractionLengthLeft = -1.0;
471 mfpKinEnergy = DBL_MAX; 465 mfpKinEnergy = DBL_MAX;
472 466
473 fParticleChange.InitializeForPostStep(track) 467 fParticleChange.InitializeForPostStep(track);
474 468
475 // Do not make anything if particle is stopp 469 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt! 470 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) 471 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478 472
479 const G4double finalT = track.GetKineticEner 473 const G4double finalT = track.GetKineticEnergy();
480 474
481 // forced process - should happen only once 475 // forced process - should happen only once per track
482 if(biasFlag) { 476 if(biasFlag) {
483 if(biasManager->ForcedInteractionRegion((G 477 if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
484 biasFlag = false; 478 biasFlag = false;
485 } 479 }
486 } 480 }
487 481
488 // check active and select model 482 // check active and select model
489 const G4double scaledEnergy = finalT*massRat 483 const G4double scaledEnergy = finalT*massRatio;
490 SelectModel(scaledEnergy, currentCoupleIndex 484 SelectModel(scaledEnergy, currentCoupleIndex);
491 if(!currentModel->IsActive(scaledEnergy)) { 485 if(!currentModel->IsActive(scaledEnergy)) { return &fParticleChange; }
492 486
493 // Integral approach 487 // Integral approach
494 if (fXSType != fEmNoIntegral) { 488 if (fXSType != fEmNoIntegral) {
495 const G4double logFinalT = 489 const G4double logFinalT =
496 track.GetDynamicParticle()->GetLogKineti 490 track.GetDynamicParticle()->GetLogKineticEnergy();
497 const G4double lx = std::max(GetCurrentLam 491 const G4double lx = std::max(GetCurrentLambda(finalT, logFinalT), 0.0);
498 #ifdef G4VERBOSE 492 #ifdef G4VERBOSE
499 if(preStepLambda < lx && 1 < verboseLevel) 493 if(preStepLambda < lx && 1 < verboseLevel) {
500 G4cout << "WARNING: for " << currentPart 494 G4cout << "WARNING: for " << currentParticle->GetParticleName()
501 << " and " << GetProcessName() << 495 << " and " << GetProcessName() << " E(MeV)= " << finalT/MeV
502 << " preLambda= " << preStepLambd 496 << " preLambda= " << preStepLambda
503 << " < " << lx << " (postLambda) 497 << " < " << lx << " (postLambda) " << G4endl;
504 } 498 }
505 #endif 499 #endif
506 // if false interaction then use new cross 500 // if false interaction then use new cross section value
507 // if both values are zero - no interactio 501 // if both values are zero - no interaction
508 if(preStepLambda*G4UniformRand() >= lx) { 502 if(preStepLambda*G4UniformRand() >= lx) {
509 return &fParticleChange; 503 return &fParticleChange;
510 } 504 }
511 } 505 }
512 506
513 // define new weight for primary and seconda 507 // define new weight for primary and secondaries
514 G4double weight = fParticleChange.GetParentW 508 G4double weight = fParticleChange.GetParentWeight();
515 if(weightFlag) { 509 if(weightFlag) {
516 weight /= biasFactor; 510 weight /= biasFactor;
517 fParticleChange.ProposeWeight(weight); 511 fParticleChange.ProposeWeight(weight);
518 } 512 }
519 513
520 #ifdef G4VERBOSE 514 #ifdef G4VERBOSE
521 if(1 < verboseLevel) { 515 if(1 < verboseLevel) {
522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 516 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV 517 << finalT/MeV
524 << " MeV; model= (" << currentModel 518 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnerg 519 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl; 520 << G4endl;
527 } 521 }
528 #endif 522 #endif
529 523
530 // sample secondaries 524 // sample secondaries
531 secParticles.clear(); 525 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticle 526 currentModel->SampleSecondaries(&secParticles,
533 currentCoupl 527 currentCouple,
534 track.GetDyn 528 track.GetDynamicParticle(),
535 (*theCuts)[c 529 (*theCuts)[currentCoupleIndex]);
536 530
537 G4int num0 = (G4int)secParticles.size(); 531 G4int num0 = (G4int)secParticles.size();
538 532
539 // splitting or Russian roulette 533 // splitting or Russian roulette
540 if(biasManager) { 534 if(biasManager) {
541 if(biasManager->SecondaryBiasingRegion((G4 535 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
542 G4double eloss = 0.0; 536 G4double eloss = 0.0;
543 weight *= biasManager->ApplySecondaryBia 537 weight *= biasManager->ApplySecondaryBiasing(
544 secParticles, track, currentModel, &fP 538 secParticles, track, currentModel, &fParticleChange, eloss,
545 (G4int)currentCoupleIndex, (*theCuts)[ 539 (G4int)currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546 step.GetPostStepPoint()->GetSafety()); 540 step.GetPostStepPoint()->GetSafety());
547 if(eloss > 0.0) { 541 if(eloss > 0.0) {
548 eloss += fParticleChange.GetLocalEnerg 542 eloss += fParticleChange.GetLocalEnergyDeposit();
549 fParticleChange.ProposeLocalEnergyDepo 543 fParticleChange.ProposeLocalEnergyDeposit(eloss);
550 } 544 }
551 } 545 }
552 } 546 }
553 547
554 // save secondaries 548 // save secondaries
555 G4int num = (G4int)secParticles.size(); 549 G4int num = (G4int)secParticles.size();
556 if(num > 0) { 550 if(num > 0) {
557 551
558 fParticleChange.SetNumberOfSecondaries(num 552 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEn 553 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime(); 554 G4double time = track.GetGlobalTime();
561 555
562 G4int n1(0), n2(0); 556 G4int n1(0), n2(0);
563 if(num0 > mainSecondaries) { << 557 if(num > mainSecondaries) {
564 currentModel->FillNumberOfSecondaries(n1 558 currentModel->FillNumberOfSecondaries(n1, n2);
565 } 559 }
566 560
567 for (G4int i=0; i<num; ++i) { 561 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i]; 562 G4DynamicParticle* dp = secParticles[i];
569 if (nullptr != dp) { 563 if (nullptr != dp) {
570 const G4ParticleDefinition* p = dp->Ge 564 const G4ParticleDefinition* p = dp->GetParticleDefinition();
571 G4double e = dp->GetKineticEnergy(); 565 G4double e = dp->GetKineticEnergy();
572 G4bool good = true; 566 G4bool good = true;
573 if(applyCuts) { 567 if(applyCuts) {
574 if (p == theGamma) { 568 if (p == theGamma) {
575 if (e < (*theCutsGamma)[currentCou 569 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576 570
577 } else if (p == theElectron) { 571 } else if (p == theElectron) {
578 if (e < (*theCutsElectron)[current 572 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579 573
580 } else if (p == thePositron) { 574 } else if (p == thePositron) {
581 if (electron_mass_c2 < (*theCutsGa 575 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582 e < (*theCutsPositron)[current 576 e < (*theCutsPositron)[currentCoupleIndex]) {
583 good = false; 577 good = false;
584 e += 2.0*electron_mass_c2; 578 e += 2.0*electron_mass_c2;
585 } 579 }
586 } 580 }
587 // added secondary if it is good 581 // added secondary if it is good
588 } 582 }
589 if (good) { 583 if (good) {
590 G4Track* t = new G4Track(dp, time, t 584 G4Track* t = new G4Track(dp, time, track.GetPosition());
591 t->SetTouchableHandle(track.GetTouch 585 t->SetTouchableHandle(track.GetTouchableHandle());
592 if (biasManager) { 586 if (biasManager) {
593 t->SetWeight(weight * biasManager- 587 t->SetWeight(weight * biasManager->GetWeight(i));
594 } else { 588 } else {
595 t->SetWeight(weight); 589 t->SetWeight(weight);
596 } 590 }
597 pParticleChange->AddSecondary(t); 591 pParticleChange->AddSecondary(t);
598 592
599 // define type of secondary 593 // define type of secondary
600 if(i < mainSecondaries) { 594 if(i < mainSecondaries) {
601 t->SetCreatorModelID(secID); 595 t->SetCreatorModelID(secID);
602 if(GetProcessSubType() == fCompton 596 if(GetProcessSubType() == fComptonScattering && p == theGamma) {
603 t->SetCreatorModelID(_ComptonGam 597 t->SetCreatorModelID(_ComptonGamma);
604 } 598 }
605 } else if(i < mainSecondaries + n1) 599 } else if(i < mainSecondaries + n1) {
606 t->SetCreatorModelID(tripletID); 600 t->SetCreatorModelID(tripletID);
607 } else if(i < mainSecondaries + n1 + 601 } else if(i < mainSecondaries + n1 + n2) {
608 t->SetCreatorModelID(_IonRecoil); 602 t->SetCreatorModelID(_IonRecoil);
609 } else { 603 } else {
610 if(i < num0) { 604 if(i < num0) {
611 if(p == theGamma) { 605 if(p == theGamma) {
612 t->SetCreatorModelID(fluoID); 606 t->SetCreatorModelID(fluoID);
613 } else { 607 } else {
614 t->SetCreatorModelID(augerID); 608 t->SetCreatorModelID(augerID);
615 } 609 }
616 } else { 610 } else {
617 t->SetCreatorModelID(biasID); << 611 t->SetCreatorModelID(secID);
618 } 612 }
619 } 613 }
620 /* 614 /*
621 G4cout << "Secondary(post step) has 615 G4cout << "Secondary(post step) has weight " << t->GetWeight()
622 << ", Ekin= " << t->GetKineti 616 << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV "
623 << GetProcessName() << " fluo 617 << GetProcessName() << " fluoID= " << fluoID
624 << " augerID= " << augerID << 618 << " augerID= " << augerID <<G4endl;
625 */ 619 */
626 } else { 620 } else {
627 delete dp; 621 delete dp;
628 edep += e; 622 edep += e;
629 } 623 }
630 } 624 }
631 } 625 }
632 fParticleChange.ProposeLocalEnergyDeposit( 626 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 } 627 }
634 628
635 if(0.0 == fParticleChange.GetProposedKinetic 629 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636 fAlive == fParticleChange.GetTrackStatus( 630 fAlive == fParticleChange.GetTrackStatus()) {
637 if(particle->GetProcessManager()->GetAtRes 631 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638 { fParticleChange.ProposeTrackStatus( 632 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639 else { fParticleChange.ProposeTrackStatus( 633 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640 } 634 }
641 635
642 return &fParticleChange; 636 return &fParticleChange;
643 } 637 }
644 638
645 //....oooOO0OOooo........oooOO0OOooo........oo 639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 640
647 G4bool G4VEmProcess::StorePhysicsTable(const G 641 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G 642 const G4String& directory,
649 G4bool 643 G4bool ascii)
650 { 644 {
651 if(!isTheMaster || part != particle) { retur 645 if(!isTheMaster || part != particle) { return true; }
652 if(G4EmTableUtil::StoreTable(this, part, the 646 if(G4EmTableUtil::StoreTable(this, part, theLambdaTable,
653 directory, "Lambda", 647 directory, "Lambda",
654 verboseLevel, a 648 verboseLevel, ascii) &&
655 G4EmTableUtil::StoreTable(this, part, the 649 G4EmTableUtil::StoreTable(this, part, theLambdaTablePrim,
656 directory, "LambdaPrim", 650 directory, "LambdaPrim",
657 verboseLevel, a 651 verboseLevel, ascii)) {
658 return true; 652 return true;
659 } 653 }
660 return false; 654 return false;
661 } 655 }
662 656
663 //....oooOO0OOooo........oooOO0OOooo........oo 657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664 658
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 659 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 cons 660 const G4String& dir,
667 G4bo 661 G4bool ascii)
668 { 662 {
669 if(!isTheMaster || part != particle) { retur 663 if(!isTheMaster || part != particle) { return true; }
670 G4bool yes = true; 664 G4bool yes = true;
671 if(buildLambdaTable) { 665 if(buildLambdaTable) {
672 yes = G4EmTableUtil::RetrieveTable(this, p 666 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTable, dir,
673 "Lambda 667 "Lambda", verboseLevel,
674 ascii, 668 ascii, splineFlag);
675 } 669 }
676 if(yes && minKinEnergyPrim < maxKinEnergy) { 670 if(yes && minKinEnergyPrim < maxKinEnergy) {
677 yes = G4EmTableUtil::RetrieveTable(this, p 671 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTablePrim, dir,
678 "Lambda 672 "LambdaPrim", verboseLevel,
679 ascii, 673 ascii, splineFlag);
680 } 674 }
681 return yes; 675 return yes;
682 } 676 }
683 677
684 //....oooOO0OOooo........oooOO0OOooo........oo 678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685 679
686 G4double G4VEmProcess::GetCrossSection(G4doubl 680 G4double G4VEmProcess::GetCrossSection(G4double kinEnergy,
687 const G 681 const G4MaterialCutsCouple* couple)
688 { 682 {
689 CurrentSetup(couple, kinEnergy); 683 CurrentSetup(couple, kinEnergy);
690 return GetCurrentLambda(kinEnergy, G4Log(kin 684 return GetCurrentLambda(kinEnergy, G4Log(kinEnergy));
691 } 685 }
692 686
693 //....oooOO0OOooo........oooOO0OOooo........oo 687 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694 688
695 G4double G4VEmProcess::GetMeanFreePath(const G 689 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696 G4doubl 690 G4double,
697 G4Force 691 G4ForceCondition* condition)
698 { 692 {
699 *condition = NotForced; 693 *condition = NotForced;
700 return G4VEmProcess::MeanFreePath(track); 694 return G4VEmProcess::MeanFreePath(track);
701 } 695 }
702 696
703 //....oooOO0OOooo........oooOO0OOooo........oo 697 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 698
705 G4double 699 G4double
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou 700 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kinEnergy,
707 G4dou 701 G4double Z, G4double A, G4double cut)
708 { 702 {
709 SelectModel(kinEnergy, currentCoupleIndex); 703 SelectModel(kinEnergy, currentCoupleIndex);
710 return (currentModel) ? 704 return (currentModel) ?
711 currentModel->ComputeCrossSectionPerAtom(c 705 currentModel->ComputeCrossSectionPerAtom(currentParticle, kinEnergy,
712 Z 706 Z, A, cut) : 0.0;
713 } 707 }
714 708
715 //....oooOO0OOooo........oooOO0OOooo........oo 709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 710
717 G4PhysicsVector* 711 G4PhysicsVector*
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate 712 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple)
719 { 713 {
720 DefineMaterial(couple); 714 DefineMaterial(couple);
721 G4PhysicsVector* newv = new G4PhysicsLogVect 715 G4PhysicsVector* newv = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy,
722 716 nLambdaBins, splineFlag);
723 return newv; 717 return newv;
724 } 718 }
725 719
726 //....oooOO0OOooo........oooOO0OOooo........oo 720 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
727 721
728 const G4Element* G4VEmProcess::GetCurrentEleme 722 const G4Element* G4VEmProcess::GetCurrentElement() const
729 { 723 {
730 return (nullptr != currentModel) ? 724 return (nullptr != currentModel) ?
731 currentModel->GetCurrentElement(currentMat 725 currentModel->GetCurrentElement(currentMaterial) : nullptr;
732 } 726 }
733 727
734 //....oooOO0OOooo........oooOO0OOooo........oo 728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735 729
736 const G4Element* G4VEmProcess::GetTargetElemen 730 const G4Element* G4VEmProcess::GetTargetElement() const
737 { 731 {
738 return (nullptr != currentModel) ? 732 return (nullptr != currentModel) ?
739 currentModel->GetCurrentElement(currentMat 733 currentModel->GetCurrentElement(currentMaterial) : nullptr;
740 } 734 }
741 735
742 //....oooOO0OOooo........oooOO0OOooo........oo 736 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743 737
744 const G4Isotope* G4VEmProcess::GetTargetIsotop 738 const G4Isotope* G4VEmProcess::GetTargetIsotope() const
745 { 739 {
746 return (nullptr != currentModel) ? 740 return (nullptr != currentModel) ?
747 currentModel->GetCurrentIsotope(GetCurrent 741 currentModel->GetCurrentIsotope(GetCurrentElement()) : nullptr;
748 } 742 }
749 743
750 //....oooOO0OOooo........oooOO0OOooo........oo 744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 745
752 void G4VEmProcess::SetCrossSectionBiasingFacto 746 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 { 747 {
754 if(f > 0.0) { 748 if(f > 0.0) {
755 biasFactor = f; 749 biasFactor = f;
756 weightFlag = flag; 750 weightFlag = flag;
757 if(1 < verboseLevel) { 751 if(1 < verboseLevel) {
758 G4cout << "### SetCrossSectionBiasingFac 752 G4cout << "### SetCrossSectionBiasingFactor: for "
759 << particle->GetParticleName() 753 << particle->GetParticleName()
760 << " and process " << GetProcessN 754 << " and process " << GetProcessName()
761 << " biasFactor= " << f << " weig 755 << " biasFactor= " << f << " weightFlag= " << flag
762 << G4endl; 756 << G4endl;
763 } 757 }
764 } 758 }
765 } 759 }
766 760
767 //....oooOO0OOooo........oooOO0OOooo........oo 761 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768 762
769 void 763 void
770 G4VEmProcess::ActivateForcedInteraction(G4doub 764 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771 G4bool 765 G4bool flag)
772 { 766 {
773 if(nullptr == biasManager) { biasManager = n 767 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
774 if(1 < verboseLevel) { 768 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: 769 G4cout << "### ActivateForcedInteraction: for "
776 << particle->GetParticleName() 770 << particle->GetParticleName()
777 << " and process " << GetProcessNam 771 << " and process " << GetProcessName()
778 << " length(mm)= " << length/mm 772 << " length(mm)= " << length/mm
779 << " in G4Region <" << r 773 << " in G4Region <" << r
780 << "> weightFlag= " << flag 774 << "> weightFlag= " << flag
781 << G4endl; 775 << G4endl;
782 } 776 }
783 weightFlag = flag; 777 weightFlag = flag;
784 biasManager->ActivateForcedInteraction(lengt 778 biasManager->ActivateForcedInteraction(length, r);
785 } 779 }
786 780
787 //....oooOO0OOooo........oooOO0OOooo........oo 781 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788 782
789 void 783 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G 784 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791 G4double factor, 785 G4double factor,
792 G4double energyLimit) 786 G4double energyLimit)
793 { 787 {
794 if (0.0 <= factor) { 788 if (0.0 <= factor) {
795 789
796 // Range cut can be applied only for e- 790 // Range cut can be applied only for e-
797 if(0.0 == factor && secondaryParticle != G 791 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798 { return; } 792 { return; }
799 793
800 if(!biasManager) { biasManager = new G4EmB 794 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801 biasManager->ActivateSecondaryBiasing(regi 795 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802 if(1 < verboseLevel) { 796 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: 797 G4cout << "### ActivateSecondaryBiasing: for "
804 << " process " << GetProcessName() 798 << " process " << GetProcessName()
805 << " factor= " << factor 799 << " factor= " << factor
806 << " in G4Region <" << region 800 << " in G4Region <" << region
807 << "> energyLimit(MeV)= " << energyLimi 801 << "> energyLimit(MeV)= " << energyLimit/MeV
808 << G4endl; 802 << G4endl;
809 } 803 }
810 } 804 }
811 } 805 }
812 806
813 //....oooOO0OOooo........oooOO0OOooo........oo 807 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 808
815 void G4VEmProcess::SetLambdaBinning(G4int n) 809 void G4VEmProcess::SetLambdaBinning(G4int n)
816 { 810 {
817 if(5 < n && n < 10000000) { 811 if(5 < n && n < 10000000) {
818 nLambdaBins = n; 812 nLambdaBins = n;
819 actBinning = true; 813 actBinning = true;
820 } else { 814 } else {
821 G4double e = (G4double)n; 815 G4double e = (G4double)n;
822 PrintWarning("SetLambdaBinning", e); 816 PrintWarning("SetLambdaBinning", e);
823 } 817 }
824 } 818 }
825 819
826 //....oooOO0OOooo........oooOO0OOooo........oo 820 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827 821
828 void G4VEmProcess::SetMinKinEnergy(G4double e) 822 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 { 823 {
830 if(1.e-3*eV < e && e < maxKinEnergy) { 824 if(1.e-3*eV < e && e < maxKinEnergy) {
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma 825 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e)
832 /G4Log(maxKinEnergy/ 826 /G4Log(maxKinEnergy/minKinEnergy));
833 minKinEnergy = e; 827 minKinEnergy = e;
834 actMinKinEnergy = true; 828 actMinKinEnergy = true;
835 } else { PrintWarning("SetMinKinEnergy", e); 829 } else { PrintWarning("SetMinKinEnergy", e); }
836 } 830 }
837 831
838 //....oooOO0OOooo........oooOO0OOooo........oo 832 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839 833
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) 834 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 { 835 {
842 if(minKinEnergy < e && e < 1.e+6*TeV) { 836 if(minKinEnergy < e && e < 1.e+6*TeV) {
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ 837 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy)
844 /G4Log(maxKinEnergy/ 838 /G4Log(maxKinEnergy/minKinEnergy));
845 maxKinEnergy = e; 839 maxKinEnergy = e;
846 actMaxKinEnergy = true; 840 actMaxKinEnergy = true;
847 } else { PrintWarning("SetMaxKinEnergy", e); 841 } else { PrintWarning("SetMaxKinEnergy", e); }
848 } 842 }
849 843
850 //....oooOO0OOooo........oooOO0OOooo........oo 844 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851 845
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl 846 void G4VEmProcess::SetMinKinEnergyPrim(G4double e)
853 { 847 {
854 if(theParameters->MinKinEnergy() <= e && 848 if(theParameters->MinKinEnergy() <= e &&
855 e <= theParameters->MaxKinEnergy()) { min 849 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; }
856 else { PrintWarning("SetMinKinEnergyPrim", e 850 else { PrintWarning("SetMinKinEnergyPrim", e); }
857 } 851 }
858 852
859 //....oooOO0OOooo........oooOO0OOooo........oo 853 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
860 854
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const 855 G4VEmProcess* G4VEmProcess::GetEmProcess(const G4String& nam)
862 { 856 {
863 return (nam == GetProcessName()) ? this : nu 857 return (nam == GetProcessName()) ? this : nullptr;
864 } 858 }
865 859
866 //....oooOO0OOooo........oooOO0OOooo........oo 860 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
867 861
868 G4double G4VEmProcess::PolarAngleLimit() const 862 G4double G4VEmProcess::PolarAngleLimit() const
869 { 863 {
870 return theParameters->MscThetaLimit(); 864 return theParameters->MscThetaLimit();
871 } 865 }
872 866
873 //....oooOO0OOooo........oooOO0OOooo........oo 867 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
874 868
875 void G4VEmProcess::PrintWarning(G4String tit, 869 void G4VEmProcess::PrintWarning(G4String tit, G4double val)
876 { 870 {
877 G4String ss = "G4VEmProcess::" + tit; 871 G4String ss = "G4VEmProcess::" + tit;
878 G4ExceptionDescription ed; 872 G4ExceptionDescription ed;
879 ed << "Parameter is out of range: " << val 873 ed << "Parameter is out of range: " << val
880 << " it will have no effect!\n" << " Pro 874 << " it will have no effect!\n" << " Process "
881 << GetProcessName() << " nbins= " << the 875 << GetProcessName() << " nbins= " << theParameters->NumberOfBins()
882 << " Emin(keV)= " << theParameters->MinKi 876 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV
883 << " Emax(GeV)= " << theParameters->MaxKi 877 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV;
884 G4Exception(ss, "em0044", JustWarning, ed); 878 G4Exception(ss, "em0044", JustWarning, ed);
885 } 879 }
886 880
887 //....oooOO0OOooo........oooOO0OOooo........oo 881 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 882
889 void G4VEmProcess::ProcessDescription(std::ost 883 void G4VEmProcess::ProcessDescription(std::ostream& out) const
890 { 884 {
891 if(nullptr != particle) { 885 if(nullptr != particle) {
892 StreamInfo(out, *particle, true); 886 StreamInfo(out, *particle, true);
893 } 887 }
894 } 888 }
895 889
896 //....oooOO0OOooo........oooOO0OOooo........oo 890 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 891