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