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