| 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 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4EmCalculator 32 // File name: G4EmCalculator
33 // 33 //
34 // Author: Vladimir Ivanchenko 34 // Author: Vladimir Ivanchenko
35 // 35 //
36 // Creation date: 28.06.2004 36 // Creation date: 28.06.2004
37 // 37 //
38 // 38 //
39 // Class Description: V.Ivanchenko & M.Novak 39 // Class Description: V.Ivanchenko & M.Novak
40 // 40 //
41 // ------------------------------------------- 41 // -------------------------------------------------------------------
42 // 42 //
43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 45
46 #include "G4EmCalculator.hh" 46 #include "G4EmCalculator.hh"
47 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
48 #include "G4LossTableManager.hh" 48 #include "G4LossTableManager.hh"
49 #include "G4EmParameters.hh" 49 #include "G4EmParameters.hh"
50 #include "G4NistManager.hh" 50 #include "G4NistManager.hh"
51 #include "G4DynamicParticle.hh" 51 #include "G4DynamicParticle.hh"
52 #include "G4VEmProcess.hh" 52 #include "G4VEmProcess.hh"
53 #include "G4VEnergyLossProcess.hh" 53 #include "G4VEnergyLossProcess.hh"
54 #include "G4VMultipleScattering.hh" 54 #include "G4VMultipleScattering.hh"
55 #include "G4Material.hh" 55 #include "G4Material.hh"
56 #include "G4MaterialCutsCouple.hh" 56 #include "G4MaterialCutsCouple.hh"
57 #include "G4ParticleDefinition.hh" 57 #include "G4ParticleDefinition.hh"
58 #include "G4ParticleTable.hh" 58 #include "G4ParticleTable.hh"
59 #include "G4IonTable.hh" 59 #include "G4IonTable.hh"
60 #include "G4PhysicsTable.hh" 60 #include "G4PhysicsTable.hh"
61 #include "G4ProductionCutsTable.hh" 61 #include "G4ProductionCutsTable.hh"
62 #include "G4ProcessManager.hh" 62 #include "G4ProcessManager.hh"
63 #include "G4ionEffectiveCharge.hh" 63 #include "G4ionEffectiveCharge.hh"
64 #include "G4RegionStore.hh" 64 #include "G4RegionStore.hh"
65 #include "G4Element.hh" 65 #include "G4Element.hh"
66 #include "G4EmCorrections.hh" 66 #include "G4EmCorrections.hh"
67 #include "G4GenericIon.hh" 67 #include "G4GenericIon.hh"
68 #include "G4ProcessVector.hh" 68 #include "G4ProcessVector.hh"
69 #include "G4Gamma.hh" 69 #include "G4Gamma.hh"
70 #include "G4Electron.hh" 70 #include "G4Electron.hh"
71 #include "G4Positron.hh" 71 #include "G4Positron.hh"
72 #include "G4EmUtility.hh" 72 #include "G4EmUtility.hh"
73 73
74 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 75
76 G4EmCalculator::G4EmCalculator() 76 G4EmCalculator::G4EmCalculator()
77 { 77 {
78 manager = G4LossTableManager::Instance(); 78 manager = G4LossTableManager::Instance();
79 nist = G4NistManager::Instance(); 79 nist = G4NistManager::Instance();
80 theParameters = G4EmParameters::Instance(); 80 theParameters = G4EmParameters::Instance();
81 corr = manager->EmCorrections(); 81 corr = manager->EmCorrections();
82 cutenergy[0] = cutenergy[1] = cutenergy[2] = 82 cutenergy[0] = cutenergy[1] = cutenergy[2] = DBL_MAX;
83 theGenericIon = G4GenericIon::GenericIon(); 83 theGenericIon = G4GenericIon::GenericIon();
84 ionEffCharge = new G4ionEffectiveCharge(); 84 ionEffCharge = new G4ionEffectiveCharge();
85 dynParticle = new G4DynamicParticle(); 85 dynParticle = new G4DynamicParticle();
86 ionTable = G4ParticleTable::GetParticle 86 ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
87 } 87 }
88 88
89 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 90
91 G4EmCalculator::~G4EmCalculator() 91 G4EmCalculator::~G4EmCalculator()
92 { 92 {
93 delete ionEffCharge; 93 delete ionEffCharge;
94 delete dynParticle; 94 delete dynParticle;
95 for (G4int i=0; i<nLocalMaterials; ++i) { 95 for (G4int i=0; i<nLocalMaterials; ++i) {
96 delete localCouples[i]; 96 delete localCouples[i];
97 } 97 }
98 } 98 }
99 99
100 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 101
102 G4double G4EmCalculator::GetDEDX(G4double kinE 102 G4double G4EmCalculator::GetDEDX(G4double kinEnergy,
103 const G4Parti 103 const G4ParticleDefinition* p,
104 const G4Mater 104 const G4Material* mat,
105 const G4Regio 105 const G4Region* region)
106 { 106 {
107 G4double res = 0.0; 107 G4double res = 0.0;
108 const G4MaterialCutsCouple* couple = FindCou 108 const G4MaterialCutsCouple* couple = FindCouple(mat, region);
109 if(nullptr != couple && UpdateParticle(p, ki 109 if(nullptr != couple && UpdateParticle(p, kinEnergy) ) {
110 res = manager->GetDEDX(p, kinEnergy, coupl 110 res = manager->GetDEDX(p, kinEnergy, couple);
111 111
112 if(isIon) { 112 if(isIon) {
113 if(FindEmModel(p, currentProcessName, ki 113 if(FindEmModel(p, currentProcessName, kinEnergy)) {
114 G4double length = CLHEP::nm; 114 G4double length = CLHEP::nm;
115 G4double eloss = res*length; 115 G4double eloss = res*length;
116 //G4cout << "### GetDEDX: E= " << kinE 116 //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res
117 // << " de= " << eloss << G4endl 117 // << " de= " << eloss << G4endl;;
118 dynParticle->SetKineticEnergy(kinEnerg 118 dynParticle->SetKineticEnergy(kinEnergy);
119 currentModel->GetChargeSquareRatio(p, 119 currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
120 currentModel->CorrectionsAlongStep(cou 120 currentModel->CorrectionsAlongStep(couple,dynParticle,length,eloss);
121 res = eloss/length; 121 res = eloss/length;
122 //G4cout << " de1= " << eloss << 122 //G4cout << " de1= " << eloss << " res1= " << res
123 // << " " << p->GetParticleName( 123 // << " " << p->GetParticleName() <<G4endl;;
124 } 124 }
125 } 125 }
126 126
127 if(verbose>0) { 127 if(verbose>0) {
128 G4cout << "G4EmCalculator::GetDEDX: E(Me 128 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
129 << " DEDX(MeV/mm)= " << res*mm/Me 129 << " DEDX(MeV/mm)= " << res*mm/MeV
130 << " DEDX(MeV*cm^2/g)= " << res*g 130 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
131 << " " << p->GetParticleName() 131 << " " << p->GetParticleName()
132 << " in " << mat->GetName() 132 << " in " << mat->GetName()
133 << " isIon= " << isIon 133 << " isIon= " << isIon
134 << G4endl; 134 << G4endl;
135 } 135 }
136 } 136 }
137 return res; 137 return res;
138 } 138 }
139 139
140 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141 141
142 G4double G4EmCalculator::GetRangeFromRestricte 142 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
143 143 const G4ParticleDefinition* p,
144 144 const G4Material* mat,
145 145 const G4Region* region)
146 { 146 {
147 G4double res = 0.0; 147 G4double res = 0.0;
148 const G4MaterialCutsCouple* couple = FindCou 148 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
149 if(couple && UpdateParticle(p, kinEnergy)) { 149 if(couple && UpdateParticle(p, kinEnergy)) {
150 res = manager->GetRangeFromRestricteDEDX(p 150 res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
151 if(verbose>1) { 151 if(verbose>1) {
152 G4cout << " G4EmCalculator::GetRangeFrom 152 G4cout << " G4EmCalculator::GetRangeFromRestrictedDEDX: E(MeV)= "
153 << kinEnergy/MeV 153 << kinEnergy/MeV
154 << " range(mm)= " << res/mm 154 << " range(mm)= " << res/mm
155 << " " << p->GetParticleName() 155 << " " << p->GetParticleName()
156 << " in " << mat->GetName() 156 << " in " << mat->GetName()
157 << G4endl; 157 << G4endl;
158 } 158 }
159 } 159 }
160 return res; 160 return res;
161 } 161 }
162 162
163 //....oooOO0OOooo........oooOO0OOooo........oo 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 164
165 G4double G4EmCalculator::GetCSDARange(G4double 165 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
166 const G4 166 const G4ParticleDefinition* p,
167 const G4 167 const G4Material* mat,
168 const G4 168 const G4Region* region)
169 { 169 {
170 G4double res = 0.0; 170 G4double res = 0.0;
171 if(!theParameters->BuildCSDARange()) { 171 if(!theParameters->BuildCSDARange()) {
172 G4ExceptionDescription ed; 172 G4ExceptionDescription ed;
173 ed << "G4EmCalculator::GetCSDARange: CSDA 173 ed << "G4EmCalculator::GetCSDARange: CSDA table is not built; "
174 << " use UI command: /process/eLoss/CSD 174 << " use UI command: /process/eLoss/CSDARange true";
175 G4Exception("G4EmCalculator::GetCSDARange" 175 G4Exception("G4EmCalculator::GetCSDARange", "em0077",
176 JustWarning, ed); 176 JustWarning, ed);
177 return res; 177 return res;
178 } 178 }
179 179
180 const G4MaterialCutsCouple* couple = FindCou 180 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
181 if(nullptr != couple && UpdateParticle(p, ki 181 if(nullptr != couple && UpdateParticle(p, kinEnergy)) {
182 res = manager->GetCSDARange(p, kinEnergy, 182 res = manager->GetCSDARange(p, kinEnergy, couple);
183 if(verbose>1) { 183 if(verbose>1) {
184 G4cout << " G4EmCalculator::GetCSDARange 184 G4cout << " G4EmCalculator::GetCSDARange: E(MeV)= " << kinEnergy/MeV
185 << " range(mm)= " << res/mm 185 << " range(mm)= " << res/mm
186 << " " << p->GetParticleName() 186 << " " << p->GetParticleName()
187 << " in " << mat->GetName() 187 << " in " << mat->GetName()
188 << G4endl; 188 << G4endl;
189 } 189 }
190 } 190 }
191 return res; 191 return res;
192 } 192 }
193 193
194 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195 195
196 G4double G4EmCalculator::GetRange(G4double kin 196 G4double G4EmCalculator::GetRange(G4double kinEnergy,
197 const G4Part 197 const G4ParticleDefinition* p,
198 const G4Mate 198 const G4Material* mat,
199 const G4Regi 199 const G4Region* region)
200 { 200 {
201 G4double res = 0.0; 201 G4double res = 0.0;
202 if(theParameters->BuildCSDARange()) { 202 if(theParameters->BuildCSDARange()) {
203 res = GetCSDARange(kinEnergy, p, mat, regi 203 res = GetCSDARange(kinEnergy, p, mat, region);
204 } else { 204 } else {
205 res = GetRangeFromRestricteDEDX(kinEnergy, 205 res = GetRangeFromRestricteDEDX(kinEnergy, p, mat, region);
206 } 206 }
207 return res; 207 return res;
208 } 208 }
209 209
210 //....oooOO0OOooo........oooOO0OOooo........oo 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 211
212 G4double G4EmCalculator::GetKinEnergy(G4double 212 G4double G4EmCalculator::GetKinEnergy(G4double range,
213 const G4 213 const G4ParticleDefinition* p,
214 const G4 214 const G4Material* mat,
215 const G4 215 const G4Region* region)
216 { 216 {
217 G4double res = 0.0; 217 G4double res = 0.0;
218 const G4MaterialCutsCouple* couple = FindCou 218 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
219 if(nullptr != couple && UpdateParticle(p, 1. 219 if(nullptr != couple && UpdateParticle(p, 1.0*GeV)) {
220 res = manager->GetEnergy(p, range, couple) 220 res = manager->GetEnergy(p, range, couple);
221 if(verbose>0) { 221 if(verbose>0) {
222 G4cout << "G4EmCalculator::GetKinEnergy: 222 G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
223 << " KinE(MeV)= " << res/MeV 223 << " KinE(MeV)= " << res/MeV
224 << " " << p->GetParticleName() 224 << " " << p->GetParticleName()
225 << " in " << mat->GetName() 225 << " in " << mat->GetName()
226 << G4endl; 226 << G4endl;
227 } 227 }
228 } 228 }
229 return res; 229 return res;
230 } 230 }
231 231
232 //....oooOO0OOooo........oooOO0OOooo........oo 232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 233
234 G4double G4EmCalculator::GetCrossSectionPerVol 234 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
235 co 235 const G4ParticleDefinition* p,
236 co 236 const G4String& processName,
237 co 237 const G4Material* mat,
238 co 238 const G4Region* region)
239 { 239 {
240 G4double res = 0.0; 240 G4double res = 0.0;
241 const G4MaterialCutsCouple* couple = FindCou 241 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
242 242
243 if(nullptr != couple && UpdateParticle(p, ki 243 if(nullptr != couple && UpdateParticle(p, kinEnergy)) {
244 if(FindEmModel(p, processName, kinEnergy)) 244 if(FindEmModel(p, processName, kinEnergy)) {
245 G4int idx = couple->GetIndex(); 245 G4int idx = couple->GetIndex();
246 G4int procType = -1; 246 G4int procType = -1;
247 FindLambdaTable(p, processName, kinEnerg 247 FindLambdaTable(p, processName, kinEnergy, procType);
248 248
249 G4VEmProcess* emproc = FindDiscreteProce 249 G4VEmProcess* emproc = FindDiscreteProcess(p, processName);
250 if(nullptr != emproc) { 250 if(nullptr != emproc) {
251 res = emproc->GetCrossSection(kinEnergy, cou 251 res = emproc->GetCrossSection(kinEnergy, couple);
252 } else if(currentLambda) { 252 } else if(currentLambda) {
253 // special tables are built for Msc mo 253 // special tables are built for Msc models
254 // procType is set in FindLambdaTable 254 // procType is set in FindLambdaTable
255 if(procType==2) { 255 if(procType==2) {
256 auto mscM = static_cast<G4VMscModel* 256 auto mscM = static_cast<G4VMscModel*>(currentModel);
257 mscM->SetCurrentCouple(couple); 257 mscM->SetCurrentCouple(couple);
258 G4double tr1Mfp = mscM->GetTransport 258 G4double tr1Mfp = mscM->GetTransportMeanFreePath(p, kinEnergy);
259 if (tr1Mfp<DBL_MAX) { 259 if (tr1Mfp<DBL_MAX) {
260 res = 1./tr1Mfp; 260 res = 1./tr1Mfp;
261 } 261 }
262 } else { 262 } else {
263 G4double e = kinEnergy*massRatio; 263 G4double e = kinEnergy*massRatio;
264 res = (((*currentLambda)[idx])->Valu 264 res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
265 } 265 }
266 } else { 266 } else {
267 res = ComputeCrossSectionPerVolume(kin 267 res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, kinEnergy);
268 } 268 }
269 if(verbose>0) { 269 if(verbose>0) {
270 G4cout << "G4EmCalculator::GetXSPerVol 270 G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
271 << " cross(cm-1)= " << res*cm 271 << " cross(cm-1)= " << res*cm
272 << " " << p->GetParticleName( 272 << " " << p->GetParticleName()
273 << " in " << mat->GetName(); 273 << " in " << mat->GetName();
274 if(verbose>1) 274 if(verbose>1)
275 G4cout << " idx= " << idx << " Esc 275 G4cout << " idx= " << idx << " Escaled((MeV)= "
276 << kinEnergy*massRatio 276 << kinEnergy*massRatio
277 << " q2= " << chargeSquare; 277 << " q2= " << chargeSquare;
278 G4cout << G4endl; 278 G4cout << G4endl;
279 } 279 }
280 } 280 }
281 } 281 }
282 return res; 282 return res;
283 } 283 }
284 284
285 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 286
287 G4double G4EmCalculator::GetShellIonisationCro 287 G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom(
288 const 288 const G4String& particle,
289 G4int 289 G4int Z,
290 G4Ato 290 G4AtomicShellEnumerator shell,
291 G4dou 291 G4double kinEnergy)
292 { 292 {
293 G4double res = 0.0; 293 G4double res = 0.0;
294 const G4ParticleDefinition* p = FindParticle 294 const G4ParticleDefinition* p = FindParticle(particle);
295 G4VAtomDeexcitation* ad = manager->AtomDeexc 295 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
296 if(nullptr != p && nullptr != ad) { 296 if(nullptr != p && nullptr != ad) {
297 res = ad->GetShellIonisationCrossSectionPe 297 res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy);
298 } 298 }
299 return res; 299 return res;
300 } 300 }
301 301
302 //....oooOO0OOooo........oooOO0OOooo........oo 302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 303
304 G4double G4EmCalculator::GetMeanFreePath(G4dou 304 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
305 const 305 const G4ParticleDefinition* p,
306 const 306 const G4String& processName,
307 const 307 const G4Material* mat,
308 const 308 const G4Region* region)
309 { 309 {
310 G4double res = DBL_MAX; 310 G4double res = DBL_MAX;
311 G4double x = GetCrossSectionPerVolume(kinEne 311 G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
312 if(x > 0.0) { res = 1.0/x; } 312 if(x > 0.0) { res = 1.0/x; }
313 if(verbose>1) { 313 if(verbose>1) {
314 G4cout << "G4EmCalculator::GetMeanFreePath 314 G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
315 << " MFP(mm)= " << res/mm 315 << " MFP(mm)= " << res/mm
316 << " " << p->GetParticleName() 316 << " " << p->GetParticleName()
317 << " in " << mat->GetName() 317 << " in " << mat->GetName()
318 << G4endl; 318 << G4endl;
319 } 319 }
320 return res; 320 return res;
321 } 321 }
322 322
323 //....oooOO0OOooo........oooOO0OOooo........oo 323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 324
325 void G4EmCalculator::PrintDEDXTable(const G4Pa 325 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
326 { 326 {
327 const G4VEnergyLossProcess* elp = manager->G 327 const G4VEnergyLossProcess* elp = manager->GetEnergyLossProcess(p);
328 G4cout << "##### DEDX Table for " << p->GetP 328 G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
329 if(nullptr != elp) G4cout << *(elp->DEDXTabl 329 if(nullptr != elp) G4cout << *(elp->DEDXTable()) << G4endl;
330 } 330 }
331 331
332 //....oooOO0OOooo........oooOO0OOooo........oo 332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333 333
334 void G4EmCalculator::PrintRangeTable(const G4P 334 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
335 { 335 {
336 const G4VEnergyLossProcess* elp = manager->G 336 const G4VEnergyLossProcess* elp = manager->GetEnergyLossProcess(p);
337 G4cout << "##### Range Table for " << p->Get 337 G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
338 if(nullptr != elp) G4cout << *(elp->RangeTab 338 if(nullptr != elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
339 } 339 }
340 340
341 //....oooOO0OOooo........oooOO0OOooo........oo 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342 342
343 void G4EmCalculator::PrintInverseRangeTable(co 343 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
344 { 344 {
345 const G4VEnergyLossProcess* elp = manager->G 345 const G4VEnergyLossProcess* elp = manager->GetEnergyLossProcess(p);
346 G4cout << "### G4EmCalculator: Inverse Range 346 G4cout << "### G4EmCalculator: Inverse Range Table for "
347 << p->GetParticleName() << G4endl; 347 << p->GetParticleName() << G4endl;
348 if(nullptr != elp) G4cout << *(elp->InverseR 348 if(nullptr != elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
349 } 349 }
350 350
351 //....oooOO0OOooo........oooOO0OOooo........oo 351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 352
353 G4double G4EmCalculator::ComputeDEDX(G4double 353 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
354 const G4P 354 const G4ParticleDefinition* p,
355 const G4S 355 const G4String& processName,
356 const G4M 356 const G4Material* mat,
357 G4d 357 G4double cut)
358 { 358 {
359 SetupMaterial(mat); 359 SetupMaterial(mat);
360 G4double res = 0.0; 360 G4double res = 0.0;
361 if(verbose > 1) { 361 if(verbose > 1) {
362 G4cout << "### G4EmCalculator::ComputeDEDX 362 G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
363 << " in " << currentMaterialName 363 << " in " << currentMaterialName
364 << " e(MeV)= " << kinEnergy/MeV << 364 << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
365 << G4endl; 365 << G4endl;
366 } 366 }
367 if(UpdateParticle(p, kinEnergy)) { 367 if(UpdateParticle(p, kinEnergy)) {
368 if(FindEmModel(p, processName, kinEnergy)) 368 if(FindEmModel(p, processName, kinEnergy)) {
369 G4double escaled = kinEnergy*massRatio; 369 G4double escaled = kinEnergy*massRatio;
370 if(nullptr != baseParticle) { 370 if(nullptr != baseParticle) {
371 res = currentModel->ComputeDEDXPerVolume(mat 371 res = currentModel->ComputeDEDXPerVolume(mat, baseParticle,
372 372 escaled, cut) * chargeSquare;
373 if(verbose > 1) { 373 if(verbose > 1) {
374 G4cout << "Particle: " << p->GetParticleNa 374 G4cout << "Particle: " << p->GetParticleName()
375 << " E(MeV)=" << kinEnergy 375 << " E(MeV)=" << kinEnergy
376 << " Base particle: " << baseParticle->Ge 376 << " Base particle: " << baseParticle->GetParticleName()
377 << " Escaled(MeV)= " << escaled 377 << " Escaled(MeV)= " << escaled
378 << " q2=" << chargeSquare << G4endl; 378 << " q2=" << chargeSquare << G4endl;
379 } 379 }
380 } else { 380 } else {
381 res = currentModel->ComputeDEDXPerVolume(mat 381 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
382 if(verbose > 1) { 382 if(verbose > 1) {
383 G4cout << "Particle: " << p->GetParticleNa 383 G4cout << "Particle: " << p->GetParticleName()
384 << " E(MeV)=" << kinEnergy << G4endl; 384 << " E(MeV)=" << kinEnergy << G4endl;
385 } 385 }
386 } 386 }
387 if(verbose > 1) { 387 if(verbose > 1) {
388 G4cout << currentModel->GetName() << ": DEDX 388 G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
389 << " DEDX(MeV*cm^2/g)= " 389 << " DEDX(MeV*cm^2/g)= "
390 << res*gram/(MeV*cm2*mat->GetDensity( 390 << res*gram/(MeV*cm2*mat->GetDensity())
391 << G4endl; 391 << G4endl;
392 } 392 }
393 // emulate smoothing procedure 393 // emulate smoothing procedure
394 if(applySmoothing && nullptr != loweMode 394 if(applySmoothing && nullptr != loweModel) {
395 G4double eth = currentModel->LowEnergyLimit( 395 G4double eth = currentModel->LowEnergyLimit();
396 G4double res0 = 0.0; 396 G4double res0 = 0.0;
397 G4double res1 = 0.0; 397 G4double res1 = 0.0;
398 if(nullptr != baseParticle) { 398 if(nullptr != baseParticle) {
399 res1 = chargeSquare* 399 res1 = chargeSquare*
400 currentModel->ComputeDEDXPerVolume(mat, 400 currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut);
401 res0 = chargeSquare* 401 res0 = chargeSquare*
402 loweModel->ComputeDEDXPerVolume(mat, bas 402 loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut);
403 } else { 403 } else {
404 res1 = currentModel->ComputeDEDXPerVolume( 404 res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
405 res0 = loweModel->ComputeDEDXPerVolume(mat 405 res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
406 } 406 }
407 if(res1 > 0.0 && escaled > 0.0) { 407 if(res1 > 0.0 && escaled > 0.0) {
408 res *= (1.0 + (res0/res1 - 1.0)*eth/escale 408 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
409 } 409 }
410 if(verbose > 1) { 410 if(verbose > 1) {
411 G4cout << "At boundary energy(MeV)= " << e 411 G4cout << "At boundary energy(MeV)= " << eth/MeV
412 << " DEDX(MeV/mm)= " << res0*mm/MeV << " 412 << " DEDX(MeV/mm)= " << res0*mm/MeV << " " << res1*mm/MeV
413 << " after correction DEDX(MeV/mm)=" << r 413 << " after correction DEDX(MeV/mm)=" << res*mm/MeV << G4endl;
414 } 414 }
415 } 415 }
416 // correction for ions 416 // correction for ions
417 if(isIon) { 417 if(isIon) {
418 const G4double length = CLHEP::nm; 418 const G4double length = CLHEP::nm;
419 if(UpdateCouple(mat, cut)) { 419 if(UpdateCouple(mat, cut)) {
420 G4double eloss = res*length; 420 G4double eloss = res*length;
421 dynParticle->SetKineticEnergy(kinEnergy); 421 dynParticle->SetKineticEnergy(kinEnergy);
422 currentModel->CorrectionsAlongStep(current 422 currentModel->CorrectionsAlongStep(currentCouple,dynParticle,
423 l 423 length,eloss);
424 res = eloss/length; 424 res = eloss/length;
425 425
426 if(verbose > 1) { 426 if(verbose > 1) {
427 G4cout << "After Corrections: DEDX(MeV/m 427 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
428 << " DEDX(MeV*cm^2/g)= " 428 << " DEDX(MeV*cm^2/g)= "
429 << res*gram/(MeV*cm2*mat->GetDensity()) 429 << res*gram/(MeV*cm2*mat->GetDensity()) << G4endl;
430 } 430 }
431 } 431 }
432 } 432 }
433 if(verbose > 0) { 433 if(verbose > 0) {
434 G4cout << "## E(MeV)= " << kinEnergy/MeV 434 G4cout << "## E(MeV)= " << kinEnergy/MeV
435 << " DEDX(MeV/mm)= " << res*mm/MeV 435 << " DEDX(MeV/mm)= " << res*mm/MeV
436 << " DEDX(MeV*cm^2/g)= " << res*gram/ 436 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
437 << " cut(MeV)= " << cut/MeV 437 << " cut(MeV)= " << cut/MeV
438 << " " << p->GetParticleName() 438 << " " << p->GetParticleName()
439 << " in " << currentMaterialName 439 << " in " << currentMaterialName
440 << " Zi^2= " << chargeSquare 440 << " Zi^2= " << chargeSquare
441 << " isIon=" << isIon 441 << " isIon=" << isIon
442 << G4endl; 442 << G4endl;
443 } 443 }
444 } 444 }
445 } 445 }
446 return res; 446 return res;
447 } 447 }
448 448
449 //....oooOO0OOooo........oooOO0OOooo........oo 449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450 450
451 G4double G4EmCalculator::ComputeElectronicDEDX 451 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
452 452 const G4ParticleDefinition* part,
453 453 const G4Material* mat,
454 454 G4double cut)
455 { 455 {
456 SetupMaterial(mat); 456 SetupMaterial(mat);
457 G4double dedx = 0.0; 457 G4double dedx = 0.0;
458 if(UpdateParticle(part, kinEnergy)) { 458 if(UpdateParticle(part, kinEnergy)) {
459 459
460 G4LossTableManager* lManager = G4LossTable 460 G4LossTableManager* lManager = G4LossTableManager::Instance();
461 const std::vector<G4VEnergyLossProcess*> v 461 const std::vector<G4VEnergyLossProcess*> vel =
462 lManager->GetEnergyLossProcessVector(); 462 lManager->GetEnergyLossProcessVector();
463 std::size_t n = vel.size(); 463 std::size_t n = vel.size();
464 464
465 //G4cout << "ComputeElectronicDEDX for " < 465 //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
466 // << " n= " << n << G4endl; 466 // << " n= " << n << G4endl;
467 467
468 for(std::size_t i=0; i<n; ++i) { 468 for(std::size_t i=0; i<n; ++i) {
469 if(vel[i]) { 469 if(vel[i]) {
470 auto p = static_cast<G4VProcess*>(vel[ 470 auto p = static_cast<G4VProcess*>(vel[i]);
471 if(ActiveForParticle(part, p)) { 471 if(ActiveForParticle(part, p)) {
472 //G4cout << "idx= " << i << " " << ( 472 //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
473 // << " " << (vel[i])->Particle()- 473 // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
474 dedx += ComputeDEDX(kinEnergy,part,( 474 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
475 } 475 }
476 } 476 }
477 } 477 }
478 } 478 }
479 return dedx; 479 return dedx;
480 } 480 }
481 481
482 //....oooOO0OOooo........oooOO0OOooo........oo 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483 483
484 G4double 484 G4double
485 G4EmCalculator::ComputeDEDXForCutInRange(G4dou 485 G4EmCalculator::ComputeDEDXForCutInRange(G4double kinEnergy,
486 const 486 const G4ParticleDefinition* part,
487 const 487 const G4Material* mat,
488 G4dou 488 G4double rangecut)
489 { 489 {
490 SetupMaterial(mat); 490 SetupMaterial(mat);
491 G4double dedx = 0.0; 491 G4double dedx = 0.0;
492 if(UpdateParticle(part, kinEnergy)) { 492 if(UpdateParticle(part, kinEnergy)) {
493 493
494 G4LossTableManager* lManager = G4LossTable 494 G4LossTableManager* lManager = G4LossTableManager::Instance();
495 const std::vector<G4VEnergyLossProcess*> v 495 const std::vector<G4VEnergyLossProcess*> vel =
496 lManager->GetEnergyLossProcessVector(); 496 lManager->GetEnergyLossProcessVector();
497 std::size_t n = vel.size(); 497 std::size_t n = vel.size();
498 498
499 if(mat != cutMaterial) { 499 if(mat != cutMaterial) {
500 cutMaterial = mat; 500 cutMaterial = mat;
501 cutenergy[0] = 501 cutenergy[0] =
502 ComputeEnergyCutFromRangeCut(rangecut, 502 ComputeEnergyCutFromRangeCut(rangecut, G4Gamma::Gamma(), mat);
503 cutenergy[1] = 503 cutenergy[1] =
504 ComputeEnergyCutFromRangeCut(rangecut, 504 ComputeEnergyCutFromRangeCut(rangecut, G4Electron::Electron(), mat);
505 cutenergy[2] = 505 cutenergy[2] =
506 ComputeEnergyCutFromRangeCut(rangecut, 506 ComputeEnergyCutFromRangeCut(rangecut, G4Positron::Positron(), mat);
507 } 507 }
508 508
509 //G4cout << "ComputeElectronicDEDX for " < 509 //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
510 // << " n= " << n << G4endl; 510 // << " n= " << n << G4endl;
511 511
512 for(std::size_t i=0; i<n; ++i) { 512 for(std::size_t i=0; i<n; ++i) {
513 if(vel[i]) { 513 if(vel[i]) {
514 auto p = static_cast<G4VProcess*>(vel[ 514 auto p = static_cast<G4VProcess*>(vel[i]);
515 if(ActiveForParticle(part, p)) { 515 if(ActiveForParticle(part, p)) {
516 //G4cout << "idx= " << i << " " << ( 516 //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
517 // << " " << (vel[i])->Particle()-> 517 // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
518 const G4ParticleDefinition* sec = (v 518 const G4ParticleDefinition* sec = (vel[i])->SecondaryParticle();
519 std::size_t idx = 0; 519 std::size_t idx = 0;
520 if(sec == G4Electron::Electron()) { 520 if(sec == G4Electron::Electron()) { idx = 1; }
521 else if(sec == G4Positron::Positron( 521 else if(sec == G4Positron::Positron()) { idx = 2; }
522 522
523 dedx += ComputeDEDX(kinEnergy,part,( 523 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),
524 mat,cutenergy[id 524 mat,cutenergy[idx]);
525 } 525 }
526 } 526 }
527 } 527 }
528 } 528 }
529 return dedx; 529 return dedx;
530 } 530 }
531 531
532 //....oooOO0OOooo........oooOO0OOooo........oo 532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
533 533
534 G4double G4EmCalculator::ComputeTotalDEDX(G4do 534 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy,
535 cons 535 const G4ParticleDefinition* part,
536 cons 536 const G4Material* mat,
537 G4do 537 G4double cut)
538 { 538 {
539 G4double dedx = ComputeElectronicDEDX(kinEne 539 G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
540 if(mass > 700.*MeV) { dedx += ComputeNuclear 540 if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
541 return dedx; 541 return dedx;
542 } 542 }
543 543
544 //....oooOO0OOooo........oooOO0OOooo........oo 544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545 545
546 G4double G4EmCalculator::ComputeNuclearDEDX(G4 546 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
547 const G4 547 const G4ParticleDefinition* p,
548 const G4 548 const G4Material* mat)
549 { 549 {
550 G4double res = 0.0; 550 G4double res = 0.0;
551 G4VEmProcess* nucst = FindDiscreteProcess(p, 551 G4VEmProcess* nucst = FindDiscreteProcess(p, "nuclearStopping");
552 if(nucst) { 552 if(nucst) {
553 G4VEmModel* mod = nucst->EmModel(); 553 G4VEmModel* mod = nucst->EmModel();
554 if(mod) { 554 if(mod) {
555 mod->SetFluctuationFlag(false); 555 mod->SetFluctuationFlag(false);
556 res = mod->ComputeDEDXPerVolume(mat, p, 556 res = mod->ComputeDEDXPerVolume(mat, p, kinEnergy);
557 } 557 }
558 } 558 }
559 559
560 if(verbose > 1) { 560 if(verbose > 1) {
561 G4cout << p->GetParticleName() << " E(MeV 561 G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
562 << " NuclearDEDX(MeV/mm)= " << res* 562 << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
563 << " NuclearDEDX(MeV*cm^2/g)= " 563 << " NuclearDEDX(MeV*cm^2/g)= "
564 << res*gram/(MeV*cm2*mat->GetDensit 564 << res*gram/(MeV*cm2*mat->GetDensity())
565 << G4endl; 565 << G4endl;
566 } 566 }
567 return res; 567 return res;
568 } 568 }
569 569
570 //....oooOO0OOooo........oooOO0OOooo........oo 570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571 571
572 G4double G4EmCalculator::ComputeCrossSectionPe 572 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
573 573 G4double kinEnergy,
574 c 574 const G4ParticleDefinition* p,
575 c 575 const G4String& processName,
576 c 576 const G4Material* mat,
577 577 G4double cut)
578 { 578 {
579 SetupMaterial(mat); 579 SetupMaterial(mat);
580 G4double res = 0.0; 580 G4double res = 0.0;
581 if(UpdateParticle(p, kinEnergy)) { 581 if(UpdateParticle(p, kinEnergy)) {
582 if(FindEmModel(p, processName, kinEnergy)) 582 if(FindEmModel(p, processName, kinEnergy)) {
583 G4double e = kinEnergy; 583 G4double e = kinEnergy;
584 G4double aCut = std::max(cut, theParamet 584 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
585 if(baseParticle) { 585 if(baseParticle) {
586 e *= kinEnergy*massRatio; 586 e *= kinEnergy*massRatio;
587 res = currentModel->CrossSectionPerVol 587 res = currentModel->CrossSectionPerVolume(
588 mat, baseParticle, e, aCut, e) * 588 mat, baseParticle, e, aCut, e) * chargeSquare;
589 } else { 589 } else {
590 res = currentModel->CrossSectionPerVol 590 res = currentModel->CrossSectionPerVolume(mat, p, e, aCut, e);
591 } 591 }
592 if(verbose>0) { 592 if(verbose>0) {
593 G4cout << "G4EmCalculator::ComputeXSPe 593 G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= "
594 << kinEnergy/MeV 594 << kinEnergy/MeV
595 << " cross(cm-1)= " << res*cm 595 << " cross(cm-1)= " << res*cm
596 << " cut(keV)= " << aCut/keV 596 << " cut(keV)= " << aCut/keV
597 << " " << p->GetParticleName( 597 << " " << p->GetParticleName()
598 << " in " << mat->GetName() 598 << " in " << mat->GetName()
599 << G4endl; 599 << G4endl;
600 } 600 }
601 } 601 }
602 } 602 }
603 return res; 603 return res;
604 } 604 }
605 605
606 //....oooOO0OOooo........oooOO0OOooo........oo 606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
607 607
608 G4double 608 G4double
609 G4EmCalculator::ComputeCrossSectionPerAtom(G4d 609 G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
610 con 610 const G4ParticleDefinition* p,
611 con 611 const G4String& processName,
612 G4d 612 G4double Z, G4double A,
613 G4d 613 G4double cut)
614 { 614 {
615 G4double res = 0.0; 615 G4double res = 0.0;
616 if(UpdateParticle(p, kinEnergy)) { 616 if(UpdateParticle(p, kinEnergy)) {
617 G4int iz = G4lrint(Z); 617 G4int iz = G4lrint(Z);
618 CheckMaterial(iz); 618 CheckMaterial(iz);
619 if(FindEmModel(p, processName, kinEnergy)) 619 if(FindEmModel(p, processName, kinEnergy)) {
620 G4double e = kinEnergy; 620 G4double e = kinEnergy;
621 G4double aCut = std::max(cut, theParamet 621 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
622 if(baseParticle) { 622 if(baseParticle) {
623 e *= kinEnergy*massRatio; 623 e *= kinEnergy*massRatio;
624 currentModel->InitialiseForElement(bas 624 currentModel->InitialiseForElement(baseParticle, iz);
625 res = currentModel->ComputeCrossSectio 625 res = currentModel->ComputeCrossSectionPerAtom(
626 baseParticle, e, Z, A, aCut) * c 626 baseParticle, e, Z, A, aCut) * chargeSquare;
627 } else { 627 } else {
628 currentModel->InitialiseForElement(p, 628 currentModel->InitialiseForElement(p, iz);
629 res = currentModel->ComputeCrossSectio 629 res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, aCut);
630 } 630 }
631 if(verbose>0) { 631 if(verbose>0) {
632 G4cout << "E(MeV)= " << kinEnergy/MeV 632 G4cout << "E(MeV)= " << kinEnergy/MeV
633 << " cross(barn)= " << res/barn 633 << " cross(barn)= " << res/barn
634 << " " << p->GetParticleName( 634 << " " << p->GetParticleName()
635 << " Z= " << Z << " A= " << A/ 635 << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
636 << " cut(keV)= " << aCut/keV 636 << " cut(keV)= " << aCut/keV
637 << G4endl; 637 << G4endl;
638 } 638 }
639 } 639 }
640 } 640 }
641 return res; 641 return res;
642 } 642 }
643 643
644 //....oooOO0OOooo........oooOO0OOooo........oo 644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645 645
646 G4double 646 G4double
647 G4EmCalculator::ComputeCrossSectionPerShell(G4 647 G4EmCalculator::ComputeCrossSectionPerShell(G4double kinEnergy,
648 co 648 const G4ParticleDefinition* p,
649 co 649 const G4String& processName,
650 G4 650 G4int Z, G4int shellIdx,
651 G4 651 G4double cut)
652 { 652 {
653 G4double res = 0.0; 653 G4double res = 0.0;
654 if(UpdateParticle(p, kinEnergy)) { 654 if(UpdateParticle(p, kinEnergy)) {
655 CheckMaterial(Z); 655 CheckMaterial(Z);
656 if(FindEmModel(p, processName, kinEnergy)) 656 if(FindEmModel(p, processName, kinEnergy)) {
657 G4double e = kinEnergy; 657 G4double e = kinEnergy;
658 G4double aCut = std::max(cut, theParamet 658 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
659 if(nullptr != baseParticle) { 659 if(nullptr != baseParticle) {
660 e *= kinEnergy*massRatio; 660 e *= kinEnergy*massRatio;
661 currentModel->InitialiseForElement(bas 661 currentModel->InitialiseForElement(baseParticle, Z);
662 res = 662 res =
663 currentModel->ComputeCrossSectionPer 663 currentModel->ComputeCrossSectionPerShell(baseParticle, Z, shellIdx,
664 664 e, aCut) * chargeSquare;
665 } else { 665 } else {
666 currentModel->InitialiseForElement(p, 666 currentModel->InitialiseForElement(p, Z);
667 res = currentModel->ComputeCrossSectio 667 res = currentModel->ComputeCrossSectionPerAtom(p, Z, shellIdx, e, aCut);
668 } 668 }
669 if(verbose>0) { 669 if(verbose>0) {
670 G4cout << "E(MeV)= " << kinEnergy/MeV 670 G4cout << "E(MeV)= " << kinEnergy/MeV
671 << " cross(barn)= " << res/barn 671 << " cross(barn)= " << res/barn
672 << " " << p->GetParticleName( 672 << " " << p->GetParticleName()
673 << " Z= " << Z << " shellIdx= 673 << " Z= " << Z << " shellIdx= " << shellIdx
674 << " cut(keV)= " << aCut/keV 674 << " cut(keV)= " << aCut/keV
675 << G4endl; 675 << G4endl;
676 } 676 }
677 } 677 }
678 } 678 }
679 return res; 679 return res;
680 } 680 }
681 681
682 //....oooOO0OOooo........oooOO0OOooo........oo 682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683 683
684 G4double 684 G4double
685 G4EmCalculator::ComputeGammaAttenuationLength( 685 G4EmCalculator::ComputeGammaAttenuationLength(G4double kinEnergy,
686 686 const G4Material* mat)
687 { 687 {
688 G4double res = 0.0; 688 G4double res = 0.0;
689 const G4ParticleDefinition* gamma = G4Gamma: 689 const G4ParticleDefinition* gamma = G4Gamma::Gamma();
690 res += ComputeCrossSectionPerVolume(kinEnerg 690 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
691 res += ComputeCrossSectionPerVolume(kinEnerg 691 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
692 res += ComputeCrossSectionPerVolume(kinEnerg 692 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
693 res += ComputeCrossSectionPerVolume(kinEnerg 693 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
694 if(res > 0.0) { res = 1.0/res; } 694 if(res > 0.0) { res = 1.0/res; }
695 return res; 695 return res;
696 } 696 }
697 697
698 //....oooOO0OOooo........oooOO0OOooo........oo 698 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
699 699
700 G4double G4EmCalculator::ComputeShellIonisatio 700 G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom(
701 const 701 const G4String& particle,
702 G4int 702 G4int Z,
703 G4Ato 703 G4AtomicShellEnumerator shell,
704 G4dou 704 G4double kinEnergy,
705 const 705 const G4Material* mat)
706 { 706 {
707 G4double res = 0.0; 707 G4double res = 0.0;
708 const G4ParticleDefinition* p = FindParticle 708 const G4ParticleDefinition* p = FindParticle(particle);
709 G4VAtomDeexcitation* ad = manager->AtomDeexc 709 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
710 if(p && ad) { 710 if(p && ad) {
711 res = ad->ComputeShellIonisationCrossSecti 711 res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell,
712 712 kinEnergy, mat);
713 } 713 }
714 return res; 714 return res;
715 } 715 }
716 716
717 //....oooOO0OOooo........oooOO0OOooo........oo 717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718 718
719 G4double G4EmCalculator::ComputeMeanFreePath(G 719 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
720 c 720 const G4ParticleDefinition* p,
721 c 721 const G4String& processName,
722 c 722 const G4Material* mat,
723 G 723 G4double cut)
724 { 724 {
725 G4double mfp = DBL_MAX; 725 G4double mfp = DBL_MAX;
726 G4double x = 726 G4double x =
727 ComputeCrossSectionPerVolume(kinEnergy, p, 727 ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
728 if(x > 0.0) { mfp = 1.0/x; } 728 if(x > 0.0) { mfp = 1.0/x; }
729 if(verbose>1) { 729 if(verbose>1) {
730 G4cout << "E(MeV)= " << kinEnergy/MeV 730 G4cout << "E(MeV)= " << kinEnergy/MeV
731 << " MFP(mm)= " << mfp/mm 731 << " MFP(mm)= " << mfp/mm
732 << " " << p->GetParticleName() 732 << " " << p->GetParticleName()
733 << " in " << mat->GetName() 733 << " in " << mat->GetName()
734 << G4endl; 734 << G4endl;
735 } 735 }
736 return mfp; 736 return mfp;
737 } 737 }
738 738
739 //....oooOO0OOooo........oooOO0OOooo........oo 739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
740 740
741 G4double G4EmCalculator::ComputeEnergyCutFromR 741 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
742 G4double range, 742 G4double range,
743 const G4ParticleDefin 743 const G4ParticleDefinition* part,
744 const G4Material* mat 744 const G4Material* mat)
745 { 745 {
746 return G4ProductionCutsTable::GetProductionC 746 return G4ProductionCutsTable::GetProductionCutsTable()->
747 ConvertRangeToEnergy(part, mat, range); 747 ConvertRangeToEnergy(part, mat, range);
748 } 748 }
749 749
750 //....oooOO0OOooo........oooOO0OOooo........oo 750 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 751
752 G4bool G4EmCalculator::UpdateParticle(const G4 752 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
753 G4double 753 G4double kinEnergy)
754 { 754 {
755 if(p != currentParticle) { 755 if(p != currentParticle) {
756 756
757 // new particle 757 // new particle
758 currentParticle = p; 758 currentParticle = p;
759 dynParticle->SetDefinition(const_cast<G4Pa 759 dynParticle->SetDefinition(const_cast<G4ParticleDefinition*>(p));
760 dynParticle->SetKineticEnergy(kinEnergy); 760 dynParticle->SetKineticEnergy(kinEnergy);
761 baseParticle = nullptr; 761 baseParticle = nullptr;
762 currentParticleName = p->GetParticleName() 762 currentParticleName = p->GetParticleName();
763 massRatio = 1.0; 763 massRatio = 1.0;
764 mass = p->GetPDGMass(); 764 mass = p->GetPDGMass();
765 chargeSquare = 1.0; 765 chargeSquare = 1.0;
766 currentProcess = manager->GetEnergyLossPr 766 currentProcess = manager->GetEnergyLossProcess(p);
767 currentProcessName = ""; 767 currentProcessName = "";
768 isIon = false; 768 isIon = false;
769 769
770 // ionisation process exist 770 // ionisation process exist
771 if(nullptr != currentProcess) { 771 if(nullptr != currentProcess) {
772 currentProcessName = currentProcess->Get 772 currentProcessName = currentProcess->GetProcessName();
773 baseParticle = currentProcess->BaseParti 773 baseParticle = currentProcess->BaseParticle();
774 if(currentProcessName == "ionIoni" && p- 774 if(currentProcessName == "ionIoni" && p->GetParticleName() != "alpha") {
775 baseParticle = theGenericIon; 775 baseParticle = theGenericIon;
776 isIon = true; 776 isIon = true;
777 } 777 }
778 778
779 // base particle is used 779 // base particle is used
780 if(nullptr != baseParticle) { 780 if(nullptr != baseParticle) {
781 massRatio = baseParticle->GetPDGMass() 781 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
782 G4double q = p->GetPDGCharge()/basePar 782 G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
783 chargeSquare = q*q; 783 chargeSquare = q*q;
784 } 784 }
785 } 785 }
786 } 786 }
787 // Effective charge for ions 787 // Effective charge for ions
788 if(isIon && nullptr != currentProcess) { 788 if(isIon && nullptr != currentProcess) {
789 chargeSquare = 789 chargeSquare =
790 corr->EffectiveChargeSquareRatio(p, curr 790 corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
791 currentProcess->SetDynamicMassCharge(massR 791 currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
792 if(verbose>1) { 792 if(verbose>1) {
793 G4cout <<"\n NewIon: massR= "<< massRati 793 G4cout <<"\n NewIon: massR= "<< massRatio << " q2= "
794 << chargeSquare << " " << currentProce 794 << chargeSquare << " " << currentProcess << G4endl;
795 } 795 }
796 } 796 }
797 return true; 797 return true;
798 } 798 }
799 799
800 //....oooOO0OOooo........oooOO0OOooo........oo 800 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
801 801
802 const G4ParticleDefinition* G4EmCalculator::Fi 802 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
803 { 803 {
804 const G4ParticleDefinition* p = nullptr; 804 const G4ParticleDefinition* p = nullptr;
805 if(name != currentParticleName) { 805 if(name != currentParticleName) {
806 p = G4ParticleTable::GetParticleTable()->F 806 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
807 if(nullptr == p) { 807 if(nullptr == p) {
808 G4cout << "### WARNING: G4EmCalculator:: 808 G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
809 << name << G4endl; 809 << name << G4endl;
810 } 810 }
811 } else { 811 } else {
812 p = currentParticle; 812 p = currentParticle;
813 } 813 }
814 return p; 814 return p;
815 } 815 }
816 816
817 //....oooOO0OOooo........oooOO0OOooo........oo 817 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
818 818
819 const G4ParticleDefinition* G4EmCalculator::Fi 819 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
820 { 820 {
821 const G4ParticleDefinition* p = ionTable->Ge 821 const G4ParticleDefinition* p = ionTable->GetIon(Z,A,0);
822 return p; 822 return p;
823 } 823 }
824 824
825 //....oooOO0OOooo........oooOO0OOooo........oo 825 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
826 826
827 const G4Material* G4EmCalculator::FindMaterial 827 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
828 { 828 {
829 if(name != currentMaterialName) { 829 if(name != currentMaterialName) {
830 SetupMaterial(G4Material::GetMaterial(name 830 SetupMaterial(G4Material::GetMaterial(name, false));
831 if(nullptr == currentMaterial) { 831 if(nullptr == currentMaterial) {
832 G4cout << "### WARNING: G4EmCalculator:: 832 G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
833 << name << G4endl; 833 << name << G4endl;
834 } 834 }
835 } 835 }
836 return currentMaterial; 836 return currentMaterial;
837 } 837 }
838 838
839 //....oooOO0OOooo........oooOO0OOooo........oo 839 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840 840
841 const G4Region* G4EmCalculator::FindRegion(con 841 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
842 { 842 {
843 return G4EmUtility::FindRegion(reg); 843 return G4EmUtility::FindRegion(reg);
844 } 844 }
845 845
846 //....oooOO0OOooo........oooOO0OOooo........oo 846 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
847 847
848 const G4MaterialCutsCouple* G4EmCalculator::Fi 848 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
849 const G4Material* 849 const G4Material* material,
850 const G4Region* re 850 const G4Region* region)
851 { 851 {
852 const G4MaterialCutsCouple* couple = nullptr 852 const G4MaterialCutsCouple* couple = nullptr;
853 SetupMaterial(material); 853 SetupMaterial(material);
854 if(nullptr != currentMaterial) { 854 if(nullptr != currentMaterial) {
855 // Access to materials 855 // Access to materials
856 const G4ProductionCutsTable* theCoupleTabl 856 const G4ProductionCutsTable* theCoupleTable=
857 G4ProductionCutsTable::GetProductionCuts 857 G4ProductionCutsTable::GetProductionCutsTable();
858 const G4Region* r = region; 858 const G4Region* r = region;
859 if(nullptr != r) { 859 if(nullptr != r) {
860 couple = theCoupleTable->GetMaterialCuts 860 couple = theCoupleTable->GetMaterialCutsCouple(material,
861 861 r->GetProductionCuts());
862 } else { 862 } else {
863 G4RegionStore* store = G4RegionStore::Ge 863 G4RegionStore* store = G4RegionStore::GetInstance();
864 std::size_t nr = store->size(); 864 std::size_t nr = store->size();
865 if(0 < nr) { 865 if(0 < nr) {
866 for(std::size_t i=0; i<nr; ++i) { 866 for(std::size_t i=0; i<nr; ++i) {
867 couple = theCoupleTable->GetMaterial 867 couple = theCoupleTable->GetMaterialCutsCouple(
868 material, ((*store)[i])->GetProduc 868 material, ((*store)[i])->GetProductionCuts());
869 if(nullptr != couple) { break; } 869 if(nullptr != couple) { break; }
870 } 870 }
871 } 871 }
872 } 872 }
873 } 873 }
874 if(nullptr == couple) { 874 if(nullptr == couple) {
875 G4ExceptionDescription ed; 875 G4ExceptionDescription ed;
876 ed << "G4EmCalculator::FindCouple: fail fo 876 ed << "G4EmCalculator::FindCouple: fail for material <"
877 << currentMaterialName << ">"; 877 << currentMaterialName << ">";
878 if(region) { ed << " and region " << regio 878 if(region) { ed << " and region " << region->GetName(); }
879 G4Exception("G4EmCalculator::FindCouple", 879 G4Exception("G4EmCalculator::FindCouple", "em0078",
880 FatalException, ed); 880 FatalException, ed);
881 } 881 }
882 return couple; 882 return couple;
883 } 883 }
884 884
885 //....oooOO0OOooo........oooOO0OOooo........oo 885 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
886 886
887 G4bool G4EmCalculator::UpdateCouple(const G4Ma 887 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
888 { 888 {
889 SetupMaterial(material); 889 SetupMaterial(material);
890 if(!currentMaterial) { return false; } 890 if(!currentMaterial) { return false; }
891 for (G4int i=0; i<nLocalMaterials; ++i) { 891 for (G4int i=0; i<nLocalMaterials; ++i) {
892 if(material == localMaterials[i] && cut == 892 if(material == localMaterials[i] && cut == localCuts[i]) {
893 currentCouple = localCouples[i]; 893 currentCouple = localCouples[i];
894 currentCoupleIndex = currentCouple->GetI 894 currentCoupleIndex = currentCouple->GetIndex();
895 currentCut = cut; 895 currentCut = cut;
896 return true; 896 return true;
897 } 897 }
898 } 898 }
899 const G4MaterialCutsCouple* cc = new G4Mater 899 const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
900 localMaterials.push_back(material); 900 localMaterials.push_back(material);
901 localCouples.push_back(cc); 901 localCouples.push_back(cc);
902 localCuts.push_back(cut); 902 localCuts.push_back(cut);
903 ++nLocalMaterials; 903 ++nLocalMaterials;
904 currentCouple = cc; 904 currentCouple = cc;
905 currentCoupleIndex = currentCouple->GetIndex 905 currentCoupleIndex = currentCouple->GetIndex();
906 currentCut = cut; 906 currentCut = cut;
907 return true; 907 return true;
908 } 908 }
909 909
910 //....oooOO0OOooo........oooOO0OOooo........oo 910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
911 911
912 void G4EmCalculator::FindLambdaTable(const G4P 912 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
913 const G4S 913 const G4String& processName,
914 G4double 914 G4double kinEnergy, G4int& proctype)
915 { 915 {
916 // Search for the process 916 // Search for the process
917 if (!currentLambda || p != lambdaParticle || 917 if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
918 lambdaName = processName; 918 lambdaName = processName;
919 currentLambda = nullptr; 919 currentLambda = nullptr;
920 lambdaParticle = p; 920 lambdaParticle = p;
921 isApplicable = false; 921 isApplicable = false;
922 922
923 const G4ParticleDefinition* part = (isIon) 923 const G4ParticleDefinition* part = (isIon) ? theGenericIon : p;
924 924
925 // Search for energy loss process 925 // Search for energy loss process
926 currentName = processName; 926 currentName = processName;
927 currentModel = nullptr; 927 currentModel = nullptr;
928 loweModel = nullptr; 928 loweModel = nullptr;
929 929
930 G4VEnergyLossProcess* elproc = FindEnLossP 930 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
931 if(nullptr != elproc) { 931 if(nullptr != elproc) {
932 currentLambda = elproc->LambdaTable(); 932 currentLambda = elproc->LambdaTable();
933 proctype = 0; 933 proctype = 0;
934 if(nullptr != currentLambda) { 934 if(nullptr != currentLambda) {
935 isApplicable = true; 935 isApplicable = true;
936 if(verbose>1) { 936 if(verbose>1) {
937 G4cout << "G4VEnergyLossProcess is f 937 G4cout << "G4VEnergyLossProcess is found out: " << currentName
938 << G4endl; 938 << G4endl;
939 } 939 }
940 } 940 }
941 curProcess = elproc; 941 curProcess = elproc;
942 return; 942 return;
943 } 943 }
944 944
945 // Search for discrete process 945 // Search for discrete process
946 G4VEmProcess* proc = FindDiscreteProcess(p 946 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
947 if(nullptr != proc) { 947 if(nullptr != proc) {
948 currentLambda = proc->LambdaTable(); 948 currentLambda = proc->LambdaTable();
949 proctype = 1; 949 proctype = 1;
950 if(nullptr != currentLambda) { 950 if(nullptr != currentLambda) {
951 isApplicable = true; 951 isApplicable = true;
952 if(verbose>1) { 952 if(verbose>1) {
953 G4cout << "G4VEmProcess is found out 953 G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
954 } 954 }
955 } 955 }
956 curProcess = proc; 956 curProcess = proc;
957 return; 957 return;
958 } 958 }
959 959
960 // Search for msc process 960 // Search for msc process
961 G4VMultipleScattering* msc = FindMscProces 961 G4VMultipleScattering* msc = FindMscProcess(part, processName);
962 if(nullptr != msc) { 962 if(nullptr != msc) {
963 currentModel = msc->SelectModel(kinEnerg 963 currentModel = msc->SelectModel(kinEnergy,0);
964 proctype = 2; 964 proctype = 2;
965 if(nullptr != currentModel) { 965 if(nullptr != currentModel) {
966 currentLambda = currentModel->GetCross 966 currentLambda = currentModel->GetCrossSectionTable();
967 if(nullptr != currentLambda) { 967 if(nullptr != currentLambda) {
968 isApplicable = true; 968 isApplicable = true;
969 if(verbose>1) { 969 if(verbose>1) {
970 G4cout << "G4VMultipleScattering i 970 G4cout << "G4VMultipleScattering is found out: " << currentName
971 << G4endl; 971 << G4endl;
972 } 972 }
973 } 973 }
974 } 974 }
975 curProcess = msc; 975 curProcess = msc;
976 } 976 }
977 } 977 }
978 } 978 }
979 979
980 //....oooOO0OOooo........oooOO0OOooo........oo 980 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
981 981
982 G4bool G4EmCalculator::FindEmModel(const G4Par 982 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
983 const G4Str 983 const G4String& processName,
984 G4 984 G4double kinEnergy)
985 { 985 {
986 isApplicable = false; 986 isApplicable = false;
987 if(nullptr == p || nullptr == currentMateria 987 if(nullptr == p || nullptr == currentMaterial) {
988 G4cout << "G4EmCalculator::FindEmModel WAR 988 G4cout << "G4EmCalculator::FindEmModel WARNING: no particle"
989 << " or materail defined; particle: 989 << " or materail defined; particle: " << p << G4endl;
990 return isApplicable; 990 return isApplicable;
991 } 991 }
992 G4String partname = p->GetParticleName(); 992 G4String partname = p->GetParticleName();
993 G4double scaledEnergy = kinEnergy*massRatio; 993 G4double scaledEnergy = kinEnergy*massRatio;
994 const G4ParticleDefinition* part = (isIon) ? 994 const G4ParticleDefinition* part = (isIon) ? theGenericIon : p;
995 995
996 if(verbose > 1) { 996 if(verbose > 1) {
997 G4cout << "## G4EmCalculator::FindEmModel 997 G4cout << "## G4EmCalculator::FindEmModel for " << partname
998 << " (type= " << p->GetParticleType 998 << " (type= " << p->GetParticleType()
999 << ") and " << processName << " at 999 << ") and " << processName << " at E(MeV)= " << scaledEnergy
1000 << G4endl; 1000 << G4endl;
1001 if(p != part) { G4cout << " GenericIon i 1001 if(p != part) { G4cout << " GenericIon is the base particle" << G4endl; }
1002 } 1002 }
1003 1003
1004 // Search for energy loss process 1004 // Search for energy loss process
1005 currentName = processName; 1005 currentName = processName;
1006 currentModel = nullptr; 1006 currentModel = nullptr;
1007 loweModel = nullptr; 1007 loweModel = nullptr;
1008 std::size_t idx = 0; 1008 std::size_t idx = 0;
1009 1009
1010 G4VEnergyLossProcess* elproc = FindEnLossPr 1010 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
1011 if(nullptr != elproc) { 1011 if(nullptr != elproc) {
1012 currentModel = elproc->SelectModelForMate 1012 currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
1013 currentModel->InitialiseForMaterial(part, 1013 currentModel->InitialiseForMaterial(part, currentMaterial);
1014 currentModel->SetupForMaterial(part, curr << 1014 currentModel->SetupForMaterial(part, currentMaterial, scaledEnergy);
1015 G4double eth = currentModel->LowEnergyLim 1015 G4double eth = currentModel->LowEnergyLimit();
1016 if(eth > 0.0) { 1016 if(eth > 0.0) {
1017 loweModel = elproc->SelectModelForMater 1017 loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1018 if(loweModel == currentModel) { loweMod 1018 if(loweModel == currentModel) { loweModel = nullptr; }
1019 else { 1019 else {
1020 loweModel->InitialiseForMaterial(part 1020 loweModel->InitialiseForMaterial(part, currentMaterial);
1021 loweModel->SetupForMaterial(part, cur 1021 loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1022 } 1022 }
1023 } 1023 }
1024 } 1024 }
1025 1025
1026 // Search for discrete process 1026 // Search for discrete process
1027 if(nullptr == currentModel) { 1027 if(nullptr == currentModel) {
1028 G4VEmProcess* proc = FindDiscreteProcess( 1028 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1029 if(nullptr != proc) { 1029 if(nullptr != proc) {
1030 currentModel = proc->SelectModelForMate 1030 currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
1031 currentModel->InitialiseForMaterial(par 1031 currentModel->InitialiseForMaterial(part, currentMaterial);
1032 currentModel->SetupForMaterial(part, cu 1032 currentModel->SetupForMaterial(part, currentMaterial, kinEnergy);
1033 G4double eth = currentModel->LowEnergyL 1033 G4double eth = currentModel->LowEnergyLimit();
1034 if(eth > 0.0) { 1034 if(eth > 0.0) {
1035 loweModel = proc->SelectModelForMater 1035 loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1036 if(loweModel == currentModel) { loweM 1036 if(loweModel == currentModel) { loweModel = nullptr; }
1037 else { 1037 else {
1038 loweModel->InitialiseForMaterial(pa 1038 loweModel->InitialiseForMaterial(part, currentMaterial);
1039 loweModel->SetupForMaterial(part, c 1039 loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1040 } 1040 }
1041 } 1041 }
1042 } 1042 }
1043 } 1043 }
1044 1044
1045 // Search for msc process 1045 // Search for msc process
1046 if(nullptr == currentModel) { 1046 if(nullptr == currentModel) {
1047 G4VMultipleScattering* proc = FindMscProc 1047 G4VMultipleScattering* proc = FindMscProcess(part, processName);
1048 if(nullptr != proc) { 1048 if(nullptr != proc) {
1049 currentModel = proc->SelectModel(kinEne 1049 currentModel = proc->SelectModel(kinEnergy, idx);
1050 loweModel = nullptr; 1050 loweModel = nullptr;
1051 } 1051 }
1052 } 1052 }
1053 if(nullptr != currentModel) { 1053 if(nullptr != currentModel) {
1054 if(loweModel == currentModel) { loweModel 1054 if(loweModel == currentModel) { loweModel = nullptr; }
1055 isApplicable = true; 1055 isApplicable = true;
1056 currentModel->InitialiseForMaterial(part, 1056 currentModel->InitialiseForMaterial(part, currentMaterial);
1057 if(loweModel) { 1057 if(loweModel) {
1058 loweModel->InitialiseForMaterial(part, 1058 loweModel->InitialiseForMaterial(part, currentMaterial);
1059 } 1059 }
1060 if(verbose > 1) { 1060 if(verbose > 1) {
1061 G4cout << " Model <" << currentModel- 1061 G4cout << " Model <" << currentModel->GetName()
1062 << "> Emin(MeV)= " << currentMod 1062 << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
1063 << " for " << part->GetParticleN 1063 << " for " << part->GetParticleName();
1064 if(nullptr != elproc) { 1064 if(nullptr != elproc) {
1065 G4cout << " and " << elproc->GetProce 1065 G4cout << " and " << elproc->GetProcessName() << " " << elproc
1066 << G4endl; 1066 << G4endl;
1067 } 1067 }
1068 if(nullptr != loweModel) { 1068 if(nullptr != loweModel) {
1069 G4cout << " LowEnergy model <" << low 1069 G4cout << " LowEnergy model <" << loweModel->GetName() << ">";
1070 } 1070 }
1071 G4cout << G4endl; 1071 G4cout << G4endl;
1072 } 1072 }
1073 } 1073 }
1074 return isApplicable; 1074 return isApplicable;
1075 } 1075 }
1076 1076
1077 //....oooOO0OOooo........oooOO0OOooo........o 1077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1078 1078
1079 G4VEnergyLossProcess* 1079 G4VEnergyLossProcess*
1080 G4EmCalculator::FindEnLossProcess(const G4Par 1080 G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
1081 const G4Str 1081 const G4String& processName)
1082 { 1082 {
1083 G4VEnergyLossProcess* proc = nullptr; 1083 G4VEnergyLossProcess* proc = nullptr;
1084 const std::vector<G4VEnergyLossProcess*> v 1084 const std::vector<G4VEnergyLossProcess*> v =
1085 manager->GetEnergyLossProcessVector(); 1085 manager->GetEnergyLossProcessVector();
1086 std::size_t n = v.size(); 1086 std::size_t n = v.size();
1087 for(std::size_t i=0; i<n; ++i) { 1087 for(std::size_t i=0; i<n; ++i) {
1088 if((v[i])->GetProcessName() == processNam 1088 if((v[i])->GetProcessName() == processName) {
1089 auto p = static_cast<G4VProcess*>(v[i]) 1089 auto p = static_cast<G4VProcess*>(v[i]);
1090 if(ActiveForParticle(part, p)) { 1090 if(ActiveForParticle(part, p)) {
1091 proc = v[i]; 1091 proc = v[i];
1092 break; 1092 break;
1093 } 1093 }
1094 } 1094 }
1095 } 1095 }
1096 return proc; 1096 return proc;
1097 } 1097 }
1098 1098
1099 //....oooOO0OOooo........oooOO0OOooo........o 1099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1100 1100
1101 G4VEmProcess* 1101 G4VEmProcess*
1102 G4EmCalculator::FindDiscreteProcess(const G4P 1102 G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
1103 const G4S 1103 const G4String& processName)
1104 { 1104 {
1105 G4VEmProcess* proc = nullptr; 1105 G4VEmProcess* proc = nullptr;
1106 auto v = manager->GetEmProcessVector(); 1106 auto v = manager->GetEmProcessVector();
1107 std::size_t n = v.size(); 1107 std::size_t n = v.size();
1108 for(std::size_t i=0; i<n; ++i) { 1108 for(std::size_t i=0; i<n; ++i) {
1109 const G4String& pName = v[i]->GetProcessN << 1109 auto pName = v[i]->GetProcessName();
1110 if(pName == "GammaGeneralProc") { 1110 if(pName == "GammaGeneralProc") {
1111 proc = v[i]->GetEmProcess(processName); 1111 proc = v[i]->GetEmProcess(processName);
1112 break; 1112 break;
1113 } else if(pName == processName) { 1113 } else if(pName == processName) {
1114 const auto p = static_cast<G4VProcess*> << 1114 auto p = static_cast<G4VProcess*>(v[i]);
1115 if(ActiveForParticle(part, p)) { 1115 if(ActiveForParticle(part, p)) {
1116 proc = v[i]; 1116 proc = v[i];
1117 break; 1117 break;
1118 } 1118 }
1119 } 1119 }
1120 } 1120 }
1121 return proc; 1121 return proc;
1122 } 1122 }
1123 1123
1124 //....oooOO0OOooo........oooOO0OOooo........o 1124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1125 1125
1126 G4VMultipleScattering* 1126 G4VMultipleScattering*
1127 G4EmCalculator::FindMscProcess(const G4Partic 1127 G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
1128 const G4String 1128 const G4String& processName)
1129 { 1129 {
1130 G4VMultipleScattering* proc = nullptr; 1130 G4VMultipleScattering* proc = nullptr;
1131 const std::vector<G4VMultipleScattering*> v 1131 const std::vector<G4VMultipleScattering*> v =
1132 manager->GetMultipleScatteringVector(); 1132 manager->GetMultipleScatteringVector();
1133 std::size_t n = v.size(); 1133 std::size_t n = v.size();
1134 for(std::size_t i=0; i<n; ++i) { 1134 for(std::size_t i=0; i<n; ++i) {
1135 if((v[i])->GetProcessName() == processNam 1135 if((v[i])->GetProcessName() == processName) {
1136 auto p = static_cast<G4VProcess*>(v[i]) 1136 auto p = static_cast<G4VProcess*>(v[i]);
1137 if(ActiveForParticle(part, p)) { 1137 if(ActiveForParticle(part, p)) {
1138 proc = v[i]; 1138 proc = v[i];
1139 break; 1139 break;
1140 } 1140 }
1141 } 1141 }
1142 } 1142 }
1143 return proc; 1143 return proc;
1144 } 1144 }
1145 1145
1146 //....oooOO0OOooo........oooOO0OOooo........o 1146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1147 1147
1148 G4VProcess* G4EmCalculator::FindProcess(const 1148 G4VProcess* G4EmCalculator::FindProcess(const G4ParticleDefinition* part,
1149 const 1149 const G4String& processName)
1150 { 1150 {
1151 G4VProcess* proc = nullptr; 1151 G4VProcess* proc = nullptr;
1152 const G4ProcessManager* procman = part->Get 1152 const G4ProcessManager* procman = part->GetProcessManager();
1153 G4ProcessVector* pv = procman->GetProcessLi 1153 G4ProcessVector* pv = procman->GetProcessList();
1154 G4int nproc = (G4int)pv->size(); 1154 G4int nproc = (G4int)pv->size();
1155 for(G4int i=0; i<nproc; ++i) { 1155 for(G4int i=0; i<nproc; ++i) {
1156 if(processName == (*pv)[i]->GetProcessNam 1156 if(processName == (*pv)[i]->GetProcessName()) {
1157 proc = (*pv)[i]; 1157 proc = (*pv)[i];
1158 break; 1158 break;
1159 } 1159 }
1160 } 1160 }
1161 return proc; 1161 return proc;
1162 } 1162 }
1163 1163
1164 //....oooOO0OOooo........oooOO0OOooo........o 1164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1165 1165
1166 G4bool G4EmCalculator::ActiveForParticle(cons 1166 G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
1167 G4VP 1167 G4VProcess* proc)
1168 { 1168 {
1169 G4ProcessManager* pm = part->GetProcessMana 1169 G4ProcessManager* pm = part->GetProcessManager();
1170 G4ProcessVector* pv = pm->GetProcessList(); 1170 G4ProcessVector* pv = pm->GetProcessList();
1171 G4int n = (G4int)pv->size(); 1171 G4int n = (G4int)pv->size();
1172 G4bool res = false; 1172 G4bool res = false;
1173 for(G4int i=0; i<n; ++i) { 1173 for(G4int i=0; i<n; ++i) {
1174 if((*pv)[i] == proc) { 1174 if((*pv)[i] == proc) {
1175 if(pm->GetProcessActivation(i)) { res = 1175 if(pm->GetProcessActivation(i)) { res = true; }
1176 break; 1176 break;
1177 } 1177 }
1178 } 1178 }
1179 return res; 1179 return res;
1180 } 1180 }
1181 1181
1182 //....oooOO0OOooo........oooOO0OOooo........o 1182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1183 1183
1184 void G4EmCalculator::SetupMaterial(const G4Ma 1184 void G4EmCalculator::SetupMaterial(const G4Material* mat)
1185 { 1185 {
1186 if(mat) { 1186 if(mat) {
1187 currentMaterial = mat; 1187 currentMaterial = mat;
1188 currentMaterialName = mat->GetName(); 1188 currentMaterialName = mat->GetName();
1189 } else { 1189 } else {
1190 currentMaterial = nullptr; 1190 currentMaterial = nullptr;
1191 currentMaterialName = ""; 1191 currentMaterialName = "";
1192 } 1192 }
1193 } 1193 }
1194 1194
1195 //....oooOO0OOooo........oooOO0OOooo........o 1195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1196 1196
1197 void G4EmCalculator::SetupMaterial(const G4St 1197 void G4EmCalculator::SetupMaterial(const G4String& mname)
1198 { 1198 {
1199 SetupMaterial(nist->FindOrBuildMaterial(mna 1199 SetupMaterial(nist->FindOrBuildMaterial(mname));
1200 } 1200 }
1201 1201
1202 //....oooOO0OOooo........oooOO0OOooo........o 1202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1203 1203
1204 void G4EmCalculator::CheckMaterial(G4int Z) 1204 void G4EmCalculator::CheckMaterial(G4int Z)
1205 { 1205 {
1206 G4bool isFound = false; 1206 G4bool isFound = false;
1207 if(nullptr != currentMaterial) { 1207 if(nullptr != currentMaterial) {
1208 G4int nn = (G4int)currentMaterial->GetNum 1208 G4int nn = (G4int)currentMaterial->GetNumberOfElements();
1209 for(G4int i=0; i<nn; ++i) { 1209 for(G4int i=0; i<nn; ++i) {
1210 if(Z == currentMaterial->GetElement(i)- 1210 if(Z == currentMaterial->GetElement(i)->GetZasInt()) {
1211 isFound = true; 1211 isFound = true;
1212 break; 1212 break;
1213 } 1213 }
1214 } 1214 }
1215 } 1215 }
1216 if(!isFound) { 1216 if(!isFound) {
1217 SetupMaterial(nist->FindOrBuildSimpleMate 1217 SetupMaterial(nist->FindOrBuildSimpleMaterial(Z));
1218 } 1218 }
1219 } 1219 }
1220 1220
1221 //....oooOO0OOooo........oooOO0OOooo........o 1221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1222 1222
1223 void G4EmCalculator::SetVerbose(G4int verb) 1223 void G4EmCalculator::SetVerbose(G4int verb)
1224 { 1224 {
1225 verbose = verb; 1225 verbose = verb;
1226 } 1226 }
1227 1227
1228 //....oooOO0OOooo........oooOO0OOooo........o 1228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1229 1229
1230 1230