| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Va 27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Valentin, M. Raine
28 // 28 //
29 // Based on the following publications 29 // Based on the following publications
30 // - Geant4 physics processes for microdo 30 // - Geant4 physics processes for microdosimetry simulation:
31 // very low energy electromagnetic models 31 // very low energy electromagnetic models for electrons in Si,
32 // NIM B, vol. 288, pp. 66 - 73, 2012. 32 // NIM B, vol. 288, pp. 66 - 73, 2012.
33 // 33 //
34 // 34 //
35 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36 36
37 37
38 #include "G4MicroElecElasticModel.hh" 38 #include "G4MicroElecElasticModel.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 #include "G4Exp.hh" 41 #include "G4Exp.hh"
42 42
43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 44
45 using namespace std; 45 using namespace std;
46 46
47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 48
49 G4MicroElecElasticModel::G4MicroElecElasticMod 49 G4MicroElecElasticModel::G4MicroElecElasticModel(const G4ParticleDefinition*,
50 const G4String& nam) 50 const G4String& nam)
51 :G4VEmModel(nam),isInitialised(false) 51 :G4VEmModel(nam),isInitialised(false)
52 { 52 {
53 nistSi = G4NistManager::Instance()->FindOrBu 53 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
54 54
55 killBelowEnergy = 16.7 * eV; // Minimum e- e 55 killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation
56 lowEnergyLimit = 0 * eV; 56 lowEnergyLimit = 0 * eV;
57 lowEnergyLimitOfModel = 5 * eV; // The model 57 lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV
58 highEnergyLimit = 100. * MeV; 58 highEnergyLimit = 100. * MeV;
59 SetLowEnergyLimit(lowEnergyLimit); 59 SetLowEnergyLimit(lowEnergyLimit);
60 SetHighEnergyLimit(highEnergyLimit); 60 SetHighEnergyLimit(highEnergyLimit);
61 61
62 verboseLevel= 0; 62 verboseLevel= 0;
63 // Verbosity scale: 63 // Verbosity scale:
64 // 0 = nothing 64 // 0 = nothing
65 // 1 = warning for energy non-conservation 65 // 1 = warning for energy non-conservation
66 // 2 = details of energy budget 66 // 2 = details of energy budget
67 // 3 = calculation of cross sections, file o 67 // 3 = calculation of cross sections, file openings, sampling of atoms
68 // 4 = entering in methods 68 // 4 = entering in methods
69 69
70 if( verboseLevel>0 ) 70 if( verboseLevel>0 )
71 { 71 {
72 G4cout << "MicroElec Elastic model is co 72 G4cout << "MicroElec Elastic model is constructed " << G4endl
73 << "Energy range: " 73 << "Energy range: "
74 << lowEnergyLimit / eV << " eV - " 74 << lowEnergyLimit / eV << " eV - "
75 << highEnergyLimit / MeV << " MeV" 75 << highEnergyLimit / MeV << " MeV"
76 << G4endl; 76 << G4endl;
77 } 77 }
78 fParticleChangeForGamma = 0; 78 fParticleChangeForGamma = 0;
79 } 79 }
80 80
81 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 82
83 G4MicroElecElasticModel::~G4MicroElecElasticMo 83 G4MicroElecElasticModel::~G4MicroElecElasticModel()
84 { 84 {
85 // For total cross section 85 // For total cross section
86 for (auto & pos : tableData) 86 for (auto & pos : tableData)
87 { 87 {
88 G4MicroElecCrossSectionDataSet* table = 88 G4MicroElecCrossSectionDataSet* table = pos.second;
89 delete table; 89 delete table;
90 } 90 }
91 91
92 // For final state 92 // For final state
93 eVecm.clear(); 93 eVecm.clear();
94 } 94 }
95 95
96 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 97
98 void G4MicroElecElasticModel::Initialise(const 98 void G4MicroElecElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
99 const G4DataVector& /*cuts*/) 99 const G4DataVector& /*cuts*/)
100 { 100 {
101 if (verboseLevel > 3) 101 if (verboseLevel > 3)
102 G4cout << "Calling G4MicroElecElasticModel 102 G4cout << "Calling G4MicroElecElasticModel::Initialise()" << G4endl;
103 103
104 // Energy limits 104 // Energy limits
105 if (LowEnergyLimit() < lowEnergyLimit) 105 if (LowEnergyLimit() < lowEnergyLimit)
106 { 106 {
107 G4cout << "G4MicroElecElasticModel: low 107 G4cout << "G4MicroElecElasticModel: low energy limit increased from " <<
108 LowEnergyLimit()/eV << " eV to " << lowEnerg 108 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
109 SetLowEnergyLimit(lowEnergyLimit); 109 SetLowEnergyLimit(lowEnergyLimit);
110 } 110 }
111 111
112 if (HighEnergyLimit() > highEnergyLimit) 112 if (HighEnergyLimit() > highEnergyLimit)
113 { 113 {
114 G4cout << "G4MicroElecElasticModel: high 114 G4cout << "G4MicroElecElasticModel: high energy limit decreased from " <<
115 HighEnergyLimit()/MeV << " MeV to " << 115 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
116 SetHighEnergyLimit(highEnergyLimit); 116 SetHighEnergyLimit(highEnergyLimit);
117 } 117 }
118 118
119 // Reading of data files 119 // Reading of data files
120 120
121 G4double scaleFactor = 1e-18 * cm * cm; 121 G4double scaleFactor = 1e-18 * cm * cm;
122 G4String fileElectron("microelec/sigma_elast 122 G4String fileElectron("microelec/sigma_elastic_e_Si");
123 123
124 G4ParticleDefinition* electronDef = G4Electr 124 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
125 G4String electron; 125 G4String electron;
126 126
127 // For total cross section 127 // For total cross section
128 electron = electronDef->GetParticleName(); 128 electron = electronDef->GetParticleName();
129 tableFile[electron] = fileElectron; 129 tableFile[electron] = fileElectron;
130 130
131 G4MicroElecCrossSectionDataSet* tableE = new 131 G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
132 tableE->LoadData(fileElectron); 132 tableE->LoadData(fileElectron);
133 tableData[electron] = tableE; 133 tableData[electron] = tableE;
134 134
135 // For final state 135 // For final state
136 const char* path = G4FindDataDir("G4LEDATA") 136 const char* path = G4FindDataDir("G4LEDATA");
137 137
138 if (!path) 138 if (!path)
139 { 139 {
140 G4Exception("G4MicroElecElasticModel::In 140 G4Exception("G4MicroElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
141 return; 141 return;
142 } 142 }
143 143
144 std::ostringstream eFullFileName; 144 std::ostringstream eFullFileName;
145 eFullFileName << path << "/microelec/sigmadi 145 eFullFileName << path << "/microelec/sigmadiff_cumulated_elastic_e_Si.dat";
146 std::ifstream eDiffCrossSection(eFullFileNam 146 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
147 147
148 if (!eDiffCrossSection) 148 if (!eDiffCrossSection)
149 G4Exception("G4MicroElecElasticModel::Init 149 G4Exception("G4MicroElecElasticModel::Initialise","em0003",FatalException,
150 "Missing data file: /microelec/sigmadiff_c 150 "Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat");
151 151
152 // Added clear for MT 152 // Added clear for MT
153 eTdummyVec.clear(); 153 eTdummyVec.clear();
154 eVecm.clear(); 154 eVecm.clear();
155 eDiffCrossSectionData.clear(); 155 eDiffCrossSectionData.clear();
156 156
157 // 157 //
158 eTdummyVec.push_back(0.); 158 eTdummyVec.push_back(0.);
159 159
160 while(!eDiffCrossSection.eof()) 160 while(!eDiffCrossSection.eof())
161 { 161 {
162 double tDummy; 162 double tDummy;
163 double eDummy; 163 double eDummy;
164 eDiffCrossSection>>tDummy>>eDummy; 164 eDiffCrossSection>>tDummy>>eDummy;
165 165
166 if (tDummy != eTdummyVec.back()) 166 if (tDummy != eTdummyVec.back())
167 { 167 {
168 eTdummyVec.push_back(tDummy); 168 eTdummyVec.push_back(tDummy);
169 eVecm[tDummy].push_back(0.); 169 eVecm[tDummy].push_back(0.);
170 } 170 }
171 171
172 eDiffCrossSection>>eDiffCrossSectionData 172 eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
173 173
174 if (eDummy != eVecm[tDummy].back()) eVec 174 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
175 } 175 }
176 // End final state 176 // End final state
177 177
178 if (verboseLevel > 2) 178 if (verboseLevel > 2)
179 G4cout << "Loaded cross section files for 179 G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
180 180
181 if( verboseLevel>0 ) 181 if( verboseLevel>0 )
182 { 182 {
183 G4cout << "MicroElec Elastic model is in 183 G4cout << "MicroElec Elastic model is initialized " << G4endl
184 << "Energy range: " 184 << "Energy range: "
185 << LowEnergyLimit() / eV << " eV - " 185 << LowEnergyLimit() / eV << " eV - "
186 << HighEnergyLimit() / MeV << " MeV" 186 << HighEnergyLimit() / MeV << " MeV"
187 << G4endl; 187 << G4endl;
188 } 188 }
189 189
190 if (isInitialised) { return; } 190 if (isInitialised) { return; }
191 fParticleChangeForGamma = GetParticleChangeF 191 fParticleChangeForGamma = GetParticleChangeForGamma();
192 isInitialised = true; 192 isInitialised = true;
193 } 193 }
194 194
195 //....oooOO0OOooo........oooOO0OOooo........oo 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 196
197 G4double G4MicroElecElasticModel::CrossSection 197 G4double G4MicroElecElasticModel::CrossSectionPerVolume(const G4Material* material,
198 const G4ParticleDefinition* p, 198 const G4ParticleDefinition* p,
199 G4double ekin, 199 G4double ekin,
200 G4double, 200 G4double,
201 G4double) 201 G4double)
202 { 202 {
203 if (verboseLevel > 3) 203 if (verboseLevel > 3)
204 G4cout << "Calling CrossSectionPerVolume() 204 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
205 205
206 // Calculate total cross section for model 206 // Calculate total cross section for model
207 G4double sigma=0; 207 G4double sigma=0;
208 G4double density = material->GetTotNbOfAtoms 208 G4double density = material->GetTotNbOfAtomsPerVolume();
209 209
210 if (material == nistSi || material->GetBaseM 210 if (material == nistSi || material->GetBaseMaterial() == nistSi)
211 { 211 {
212 const G4String& particleName = p->GetPar 212 const G4String& particleName = p->GetParticleName();
213 213
214 if (ekin < highEnergyLimit) 214 if (ekin < highEnergyLimit)
215 { 215 {
216 //SI : XS must not be zero otherwise sampl 216 //SI : XS must not be zero otherwise sampling of secondaries method ignored
217 if (ekin < killBelowEnergy) return DBL_MAX 217 if (ekin < killBelowEnergy) return DBL_MAX;
218 // 218 //
219 219
220 auto pos = tableData.find(particleName); 220 auto pos = tableData.find(particleName);
221 if (pos != tableData.end()) 221 if (pos != tableData.end())
222 { 222 {
223 G4MicroElecCrossSectionDataSet* table 223 G4MicroElecCrossSectionDataSet* table = pos->second;
224 if (table != nullptr) 224 if (table != nullptr)
225 { 225 {
226 sigma = table->FindValue(ekin); 226 sigma = table->FindValue(ekin);
227 } 227 }
228 } 228 }
229 else 229 else
230 { 230 {
231 G4Exception("G4MicroElecElasticModel:: 231 G4Exception("G4MicroElecElasticModel::ComputeCrossSectionPerVolume","em0002",
232 FatalException,"Model not applicable t 232 FatalException,"Model not applicable to particle type.");
233 } 233 }
234 } 234 }
235 235
236 if (verboseLevel > 3) 236 if (verboseLevel > 3)
237 { 237 {
238 G4cout << "---> Kinetic energy(eV)=" << ek 238 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
239 G4cout << " - Cross section per Si atom (c 239 G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
240 G4cout << " - Cross section per Si atom (c 240 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
241 } 241 }
242 } 242 }
243 return sigma*density; 243 return sigma*density;
244 } 244 }
245 245
246 //....oooOO0OOooo........oooOO0OOooo........oo 246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247 247
248 void G4MicroElecElasticModel::SampleSecondarie 248 void G4MicroElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
249 const G4MaterialCutsCouple* /*coup 249 const G4MaterialCutsCouple* /*couple*/,
250 const G4DynamicParticle* aDynamicE 250 const G4DynamicParticle* aDynamicElectron,
251 G4double, 251 G4double,
252 G4double) 252 G4double)
253 { 253 {
254 254
255 if (verboseLevel > 3) 255 if (verboseLevel > 3)
256 G4cout << "Calling SampleSecondaries() of 256 G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
257 257
258 G4double electronEnergy0 = aDynamicElectron- 258 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
259 259
260 if (electronEnergy0 < killBelowEnergy) 260 if (electronEnergy0 < killBelowEnergy)
261 { 261 {
262 fParticleChangeForGamma->SetProposedKine 262 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
263 fParticleChangeForGamma->ProposeTrackSta 263 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
264 fParticleChangeForGamma->ProposeLocalEne 264 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
265 return ; 265 return ;
266 } 266 }
267 267
268 if (electronEnergy0>= killBelowEnergy && ele 268 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
269 { 269 {
270 G4double cosTheta = RandomizeCosTheta(el 270 G4double cosTheta = RandomizeCosTheta(electronEnergy0);
271 G4double phi = 2. * pi * G4UniformRand() 271 G4double phi = 2. * pi * G4UniformRand();
272 G4ThreeVector zVers = aDynamicElectron-> 272 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
273 G4ThreeVector xVers = zVers.orthogonal() 273 G4ThreeVector xVers = zVers.orthogonal();
274 G4ThreeVector yVers = zVers.cross(xVers) 274 G4ThreeVector yVers = zVers.cross(xVers);
275 275
276 G4double xDir = std::sqrt(1. - cosTheta* 276 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
277 G4double yDir = xDir; 277 G4double yDir = xDir;
278 xDir *= std::cos(phi); 278 xDir *= std::cos(phi);
279 yDir *= std::sin(phi); 279 yDir *= std::sin(phi);
280 280
281 G4ThreeVector zPrimeVers((xDir*xVers + y 281 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
282 282
283 fParticleChangeForGamma->ProposeMomentum 283 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
284 fParticleChangeForGamma->SetProposedKine 284 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
285 } 285 }
286 } 286 }
287 287
288 //....oooOO0OOooo........oooOO0OOooo........oo 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 289
290 G4double G4MicroElecElasticModel::Theta 290 G4double G4MicroElecElasticModel::Theta
291 (G4ParticleDefinition * particleDefinition, G4 291 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
292 { 292 {
293 G4double theta = 0.; 293 G4double theta = 0.;
294 G4double valueT1 = 0; 294 G4double valueT1 = 0;
295 G4double valueT2 = 0; 295 G4double valueT2 = 0;
296 G4double valueE21 = 0; 296 G4double valueE21 = 0;
297 G4double valueE22 = 0; 297 G4double valueE22 = 0;
298 G4double valueE12 = 0; 298 G4double valueE12 = 0;
299 G4double valueE11 = 0; 299 G4double valueE11 = 0;
300 G4double xs11 = 0; 300 G4double xs11 = 0;
301 G4double xs12 = 0; 301 G4double xs12 = 0;
302 G4double xs21 = 0; 302 G4double xs21 = 0;
303 G4double xs22 = 0; 303 G4double xs22 = 0;
304 304
305 if (particleDefinition == G4Electron::Electr 305 if (particleDefinition == G4Electron::ElectronDefinition())
306 { 306 {
307 auto t2 = std::upper_bound(eTdummyVec.be 307 auto t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
308 auto t1 = t2-1; 308 auto t1 = t2-1;
309 auto e12 = std::upper_bound(eVecm[(*t1)] 309 auto e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
310 auto e11 = e12-1; 310 auto e11 = e12-1;
311 auto e22 = std::upper_bound(eVecm[(*t2)] 311 auto e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
312 auto e21 = e22-1; 312 auto e21 = e22-1;
313 313
314 valueT1 =*t1; 314 valueT1 =*t1;
315 valueT2 =*t2; 315 valueT2 =*t2;
316 valueE21 =*e21; 316 valueE21 =*e21;
317 valueE22 =*e22; 317 valueE22 =*e22;
318 valueE12 =*e12; 318 valueE12 =*e12;
319 valueE11 =*e11; 319 valueE11 =*e11;
320 320
321 xs11 = eDiffCrossSectionData[valueT1][va 321 xs11 = eDiffCrossSectionData[valueT1][valueE11];
322 xs12 = eDiffCrossSectionData[valueT1][va 322 xs12 = eDiffCrossSectionData[valueT1][valueE12];
323 xs21 = eDiffCrossSectionData[valueT2][va 323 xs21 = eDiffCrossSectionData[valueT2][valueE21];
324 xs22 = eDiffCrossSectionData[valueT2][va 324 xs22 = eDiffCrossSectionData[valueT2][valueE22];
325 } 325 }
326 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 326 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
327 327
328 theta = QuadInterpolator( valueE11, valueE1 328 theta = QuadInterpolator( valueE11, valueE12,
329 valueE21, valueE22, 329 valueE21, valueE22,
330 xs11, xs12, 330 xs11, xs12,
331 xs21, xs22, 331 xs21, xs22,
332 valueT1, valueT2, 332 valueT1, valueT2,
333 k, integrDiff ); 333 k, integrDiff );
334 334
335 return theta; 335 return theta;
336 } 336 }
337 337
338 //....oooOO0OOooo........oooOO0OOooo........oo 338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 339
340 G4double G4MicroElecElasticModel::LinLogInterp 340 G4double G4MicroElecElasticModel::LinLogInterpolate(G4double e1,
341 G4double e2, 341 G4double e2,
342 G4double e, 342 G4double e,
343 G4double xs1, 343 G4double xs1,
344 G4double xs2) 344 G4double xs2)
345 { 345 {
346 G4double d1 = std::log(xs1); 346 G4double d1 = std::log(xs1);
347 G4double d2 = std::log(xs2); 347 G4double d2 = std::log(xs2);
348 G4double value = G4Exp(d1 + (d2 - d1)*(e - e 348 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
349 return value; 349 return value;
350 } 350 }
351 351
352 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353 353
354 G4double G4MicroElecElasticModel::LinLinInterp 354 G4double G4MicroElecElasticModel::LinLinInterpolate(G4double e1,
355 G4double e2, 355 G4double e2,
356 G4double e, 356 G4double e,
357 G4double xs1, 357 G4double xs1,
358 G4double xs2) 358 G4double xs2)
359 { 359 {
360 G4double d1 = xs1; 360 G4double d1 = xs1;
361 G4double d2 = xs2; 361 G4double d2 = xs2;
362 G4double value = (d1 + (d2 - d1)*(e - e1)/ ( 362 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
363 return value; 363 return value;
364 } 364 }
365 365
366 //....oooOO0OOooo........oooOO0OOooo........oo 366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367 367
368 G4double G4MicroElecElasticModel::LogLogInterp 368 G4double G4MicroElecElasticModel::LogLogInterpolate(G4double e1,
369 G4double e2, 369 G4double e2,
370 G4double e, 370 G4double e,
371 G4double xs1, 371 G4double xs1,
372 G4double xs2) 372 G4double xs2)
373 { 373 {
374 G4double a = (std::log10(xs2)-std::log10(xs1 374 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
375 G4double b = std::log10(xs2) - a*std::log10( 375 G4double b = std::log10(xs2) - a*std::log10(e2);
376 G4double sigma = a*std::log10(e) + b; 376 G4double sigma = a*std::log10(e) + b;
377 G4double value = (std::pow(10.,sigma)); 377 G4double value = (std::pow(10.,sigma));
378 return value; 378 return value;
379 } 379 }
380 380
381 //....oooOO0OOooo........oooOO0OOooo........oo 381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
382 382
383 G4double G4MicroElecElasticModel::QuadInterpol 383 G4double G4MicroElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
384 G4double e21, G4double e22, 384 G4double e21, G4double e22,
385 G4double xs11, G4double xs12, 385 G4double xs11, G4double xs12,
386 G4double xs21, G4double xs22, 386 G4double xs21, G4double xs22,
387 G4double t1, G4double t2, 387 G4double t1, G4double t2,
388 G4double t, G4double e) 388 G4double t, G4double e)
389 { 389 {
390 // Log-Log 390 // Log-Log
391 /* 391 /*
392 G4double interpolatedvalue1 = LogLogInterp 392 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
393 G4double interpolatedvalue2 = LogLogInterp 393 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
394 G4double value = LogLogInterpolate(t1, t2, 394 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
395 395
396 396
397 // Lin-Log 397 // Lin-Log
398 G4double interpolatedvalue1 = LinLogInterp 398 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
399 G4double interpolatedvalue2 = LinLogInterp 399 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
400 G4double value = LinLogInterpolate(t1, t2, 400 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
401 */ 401 */
402 402
403 // Lin-Lin 403 // Lin-Lin
404 G4double interpolatedvalue1 = LinLinInterpol 404 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
405 G4double interpolatedvalue2 = LinLinInterpol 405 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
406 G4double value = LinLinInterpolate(t1, t2, t 406 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
407 407
408 return value; 408 return value;
409 } 409 }
410 410
411 //....oooOO0OOooo........oooOO0OOooo........oo 411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412 412
413 G4double G4MicroElecElasticModel::RandomizeCos 413 G4double G4MicroElecElasticModel::RandomizeCosTheta(G4double k)
414 { 414 {
415 G4double integrdiff=0; 415 G4double integrdiff=0;
416 G4double uniformRand=G4UniformRand(); 416 G4double uniformRand=G4UniformRand();
417 integrdiff = uniformRand; 417 integrdiff = uniformRand;
418 418
419 G4double theta=0.; 419 G4double theta=0.;
420 G4double cosTheta=0.; 420 G4double cosTheta=0.;
421 theta = Theta(G4Electron::ElectronDefinition 421 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
422 422
423 cosTheta= std::cos(theta*pi/180); 423 cosTheta= std::cos(theta*pi/180);
424 424
425 return cosTheta; 425 return cosTheta;
426 } 426 }
427 427