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