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