| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Va
28 //
29 // Based on the following publications
30 // - Geant4 physics processes for microdo
31 // very low energy electromagnetic models
32 // NIM B, vol. 288, pp. 66 - 73, 2012.
33 //
34 //
35 //....oooOO0OOooo........oooOO0OOooo........oo
36
37
38 #include "G4MicroElecElasticModel.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Exp.hh"
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 using namespace std;
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4MicroElecElasticModel::G4MicroElecElasticMod
50 const G4String& nam)
51 :G4VEmModel(nam),isInitialised(false)
52 {
53 nistSi = G4NistManager::Instance()->FindOrBu
54
55 killBelowEnergy = 16.7 * eV; // Minimum e- e
56 lowEnergyLimit = 0 * eV;
57 lowEnergyLimitOfModel = 5 * eV; // The model
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 o
68 // 4 = entering in methods
69
70 if( verboseLevel>0 )
71 {
72 G4cout << "MicroElec Elastic model is co
73 << "Energy range: "
74 << lowEnergyLimit / eV << " eV - "
75 << highEnergyLimit / MeV << " MeV"
76 << G4endl;
77 }
78 fParticleChangeForGamma = 0;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4MicroElecElasticModel::~G4MicroElecElasticMo
84 {
85 // For total cross section
86 for (auto & pos : tableData)
87 {
88 G4MicroElecCrossSectionDataSet* table =
89 delete table;
90 }
91
92 // For final state
93 eVecm.clear();
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97
98 void G4MicroElecElasticModel::Initialise(const
99 const G4DataVector& /*cuts*/)
100 {
101 if (verboseLevel > 3)
102 G4cout << "Calling G4MicroElecElasticModel
103
104 // Energy limits
105 if (LowEnergyLimit() < lowEnergyLimit)
106 {
107 G4cout << "G4MicroElecElasticModel: low
108 LowEnergyLimit()/eV << " eV to " << lowEnerg
109 SetLowEnergyLimit(lowEnergyLimit);
110 }
111
112 if (HighEnergyLimit() > highEnergyLimit)
113 {
114 G4cout << "G4MicroElecElasticModel: high
115 HighEnergyLimit()/MeV << " MeV to " <<
116 SetHighEnergyLimit(highEnergyLimit);
117 }
118
119 // Reading of data files
120
121 G4double scaleFactor = 1e-18 * cm * cm;
122 G4String fileElectron("microelec/sigma_elast
123
124 G4ParticleDefinition* electronDef = G4Electr
125 G4String electron;
126
127 // For total cross section
128 electron = electronDef->GetParticleName();
129 tableFile[electron] = fileElectron;
130
131 G4MicroElecCrossSectionDataSet* tableE = new
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::In
141 return;
142 }
143
144 std::ostringstream eFullFileName;
145 eFullFileName << path << "/microelec/sigmadi
146 std::ifstream eDiffCrossSection(eFullFileNam
147
148 if (!eDiffCrossSection)
149 G4Exception("G4MicroElecElasticModel::Init
150 "Missing data file: /microelec/sigmadiff_c
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
173
174 if (eDummy != eVecm[tDummy].back()) eVec
175 }
176 // End final state
177
178 if (verboseLevel > 2)
179 G4cout << "Loaded cross section files for
180
181 if( verboseLevel>0 )
182 {
183 G4cout << "MicroElec Elastic model is in
184 << "Energy range: "
185 << LowEnergyLimit() / eV << " eV - "
186 << HighEnergyLimit() / MeV << " MeV"
187 << G4endl;
188 }
189
190 if (isInitialised) { return; }
191 fParticleChangeForGamma = GetParticleChangeF
192 isInitialised = true;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oo
196
197 G4double G4MicroElecElasticModel::CrossSection
198 const G4ParticleDefinition* p,
199 G4double ekin,
200 G4double,
201 G4double)
202 {
203 if (verboseLevel > 3)
204 G4cout << "Calling CrossSectionPerVolume()
205
206 // Calculate total cross section for model
207 G4double sigma=0;
208 G4double density = material->GetTotNbOfAtoms
209
210 if (material == nistSi || material->GetBaseM
211 {
212 const G4String& particleName = p->GetPar
213
214 if (ekin < highEnergyLimit)
215 {
216 //SI : XS must not be zero otherwise sampl
217 if (ekin < killBelowEnergy) return DBL_MAX
218 //
219
220 auto pos = tableData.find(particleName);
221 if (pos != tableData.end())
222 {
223 G4MicroElecCrossSectionDataSet* table
224 if (table != nullptr)
225 {
226 sigma = table->FindValue(ekin);
227 }
228 }
229 else
230 {
231 G4Exception("G4MicroElecElasticModel::
232 FatalException,"Model not applicable t
233 }
234 }
235
236 if (verboseLevel > 3)
237 {
238 G4cout << "---> Kinetic energy(eV)=" << ek
239 G4cout << " - Cross section per Si atom (c
240 G4cout << " - Cross section per Si atom (c
241 }
242 }
243 return sigma*density;
244 }
245
246 //....oooOO0OOooo........oooOO0OOooo........oo
247
248 void G4MicroElecElasticModel::SampleSecondarie
249 const G4MaterialCutsCouple* /*coup
250 const G4DynamicParticle* aDynamicE
251 G4double,
252 G4double)
253 {
254
255 if (verboseLevel > 3)
256 G4cout << "Calling SampleSecondaries() of
257
258 G4double electronEnergy0 = aDynamicElectron-
259
260 if (electronEnergy0 < killBelowEnergy)
261 {
262 fParticleChangeForGamma->SetProposedKine
263 fParticleChangeForGamma->ProposeTrackSta
264 fParticleChangeForGamma->ProposeLocalEne
265 return ;
266 }
267
268 if (electronEnergy0>= killBelowEnergy && ele
269 {
270 G4double cosTheta = RandomizeCosTheta(el
271 G4double phi = 2. * pi * G4UniformRand()
272 G4ThreeVector zVers = aDynamicElectron->
273 G4ThreeVector xVers = zVers.orthogonal()
274 G4ThreeVector yVers = zVers.cross(xVers)
275
276 G4double xDir = std::sqrt(1. - cosTheta*
277 G4double yDir = xDir;
278 xDir *= std::cos(phi);
279 yDir *= std::sin(phi);
280
281 G4ThreeVector zPrimeVers((xDir*xVers + y
282
283 fParticleChangeForGamma->ProposeMomentum
284 fParticleChangeForGamma->SetProposedKine
285 }
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 G4double G4MicroElecElasticModel::Theta
291 (G4ParticleDefinition * particleDefinition, G4
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::Electr
306 {
307 auto t2 = std::upper_bound(eTdummyVec.be
308 auto t1 = t2-1;
309 auto e12 = std::upper_bound(eVecm[(*t1)]
310 auto e11 = e12-1;
311 auto e22 = std::upper_bound(eVecm[(*t2)]
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][va
322 xs12 = eDiffCrossSectionData[valueT1][va
323 xs21 = eDiffCrossSectionData[valueT2][va
324 xs22 = eDiffCrossSectionData[valueT2][va
325 }
326 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
327
328 theta = QuadInterpolator( valueE11, valueE1
329 valueE21, valueE22,
330 xs11, xs12,
331 xs21, xs22,
332 valueT1, valueT2,
333 k, integrDiff );
334
335 return theta;
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oo
339
340 G4double G4MicroElecElasticModel::LinLogInterp
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 - e
349 return value;
350 }
351
352 //....oooOO0OOooo........oooOO0OOooo........oo
353
354 G4double G4MicroElecElasticModel::LinLinInterp
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)/ (
363 return value;
364 }
365
366 //....oooOO0OOooo........oooOO0OOooo........oo
367
368 G4double G4MicroElecElasticModel::LogLogInterp
369 G4double e2,
370 G4double e,
371 G4double xs1,
372 G4double xs2)
373 {
374 G4double a = (std::log10(xs2)-std::log10(xs1
375 G4double b = std::log10(xs2) - a*std::log10(
376 G4double sigma = a*std::log10(e) + b;
377 G4double value = (std::pow(10.,sigma));
378 return value;
379 }
380
381 //....oooOO0OOooo........oooOO0OOooo........oo
382
383 G4double G4MicroElecElasticModel::QuadInterpol
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 = LogLogInterp
393 G4double interpolatedvalue2 = LogLogInterp
394 G4double value = LogLogInterpolate(t1, t2,
395
396
397 // Lin-Log
398 G4double interpolatedvalue1 = LinLogInterp
399 G4double interpolatedvalue2 = LinLogInterp
400 G4double value = LinLogInterpolate(t1, t2,
401 */
402
403 // Lin-Lin
404 G4double interpolatedvalue1 = LinLinInterpol
405 G4double interpolatedvalue2 = LinLinInterpol
406 G4double value = LinLinInterpolate(t1, t2, t
407
408 return value;
409 }
410
411 //....oooOO0OOooo........oooOO0OOooo........oo
412
413 G4double G4MicroElecElasticModel::RandomizeCos
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
422
423 cosTheta= std::cos(theta*pi/180);
424
425 return cosTheta;
426 }
427