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