Geant4  10.01.p02
G4MicroElecElasticModel.cc
Go to the documentation of this file.
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 
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 
44 using namespace std;
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
49  const G4String& nam)
50 :G4VEmModel(nam),isInitialised(false)
51 {
53 
54  killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation
55  lowEnergyLimit = 0 * eV;
56  lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV
57  highEnergyLimit = 100. * MeV;
60 
61  verboseLevel= 0;
62  // Verbosity scale:
63  // 0 = nothing
64  // 1 = warning for energy non-conservation
65  // 2 = details of energy budget
66  // 3 = calculation of cross sections, file openings, sampling of atoms
67  // 4 = entering in methods
68 
69  if( verboseLevel>0 )
70  {
71  G4cout << "MicroElec Elastic model is constructed " << G4endl
72  << "Energy range: "
73  << lowEnergyLimit / eV << " eV - "
74  << highEnergyLimit / MeV << " MeV"
75  << G4endl;
76  }
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
83 {
84  // For total cross section
85 
86  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
87  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
88  {
89  G4MicroElecCrossSectionDataSet* table = pos->second;
90  delete table;
91  }
92 
93  // For final state
94 
95  eVecm.clear();
96 
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
102  const G4DataVector& /*cuts*/)
103 {
104 
105  if (verboseLevel > 3)
106  G4cout << "Calling G4MicroElecElasticModel::Initialise()" << G4endl;
107 
108  // Energy limits
109 
111  {
112  G4cout << "G4MicroElecElasticModel: low energy limit increased from " <<
113  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
115  }
116 
118  {
119  G4cout << "G4MicroElecElasticModel: high energy limit decreased from " <<
120  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
122  }
123 
124  // Reading of data files
125 
126  G4double scaleFactor = 1e-18 * cm * cm;
127 
128  G4String fileElectron("microelec/sigma_elastic_e_Si");
129 
132 
133  // For total cross section
134 
135  electron = electronDef->GetParticleName();
136 
137  tableFile[electron] = fileElectron;
138 
140  tableE->LoadData(fileElectron);
141  tableData[electron] = tableE;
142 
143  // For final state
144 
145  char *path = getenv("G4LEDATA");
146 
147  if (!path)
148  {
149  G4Exception("G4MicroElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
150  return;
151  }
152 
153  std::ostringstream eFullFileName;
154  eFullFileName << path << "/microelec/sigmadiff_cumulated_elastic_e_Si.dat";
155  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
156 
157  if (!eDiffCrossSection)
158  G4Exception("G4MicroElecElasticModel::Initialise","em0003",FatalException,"Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat");
159 
160 
161  // October 21th, 2014 - Melanie Raine
162  // Added clear for MT
163 
164  eTdummyVec.clear();
165  eVecm.clear();
166  eDiffCrossSectionData.clear();
167 
168  //
169 
170 
171  eTdummyVec.push_back(0.);
172 
173  while(!eDiffCrossSection.eof())
174  {
175  double tDummy;
176  double eDummy;
177  eDiffCrossSection>>tDummy>>eDummy;
178 
179  // SI : mandatory eVecm initialization
180 
181  if (tDummy != eTdummyVec.back())
182  {
183  eTdummyVec.push_back(tDummy);
184  eVecm[tDummy].push_back(0.);
185  }
186 
187  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
188 
189  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
190 
191  }
192 
193  // End final state
194 
195  if (verboseLevel > 2)
196  G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
197 
198  if( verboseLevel>0 )
199  {
200  G4cout << "MicroElec Elastic model is initialized " << G4endl
201  << "Energy range: "
202  << LowEnergyLimit() / eV << " eV - "
203  << HighEnergyLimit() / MeV << " MeV"
204  << G4endl;
205  }
206 
207  if (isInitialised) { return; }
209  isInitialised = true;
210 
211 }
212 
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214 
216  const G4ParticleDefinition* p,
217  G4double ekin,
218  G4double,
219  G4double)
220 {
221  if (verboseLevel > 3)
222  G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
223 
224  // Calculate total cross section for model
225 
226  G4double sigma=0;
227 
229 
230  if (material == nistSi || material->GetBaseMaterial() == nistSi)
231  {
232  const G4String& particleName = p->GetParticleName();
233 
234  if (ekin < highEnergyLimit)
235  {
236  //SI : XS must not be zero otherwise sampling of secondaries method ignored
237  if (ekin < killBelowEnergy) return DBL_MAX;
238  //
239 
240  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
241  pos = tableData.find(particleName);
242 
243  if (pos != tableData.end())
244  {
245  G4MicroElecCrossSectionDataSet* table = pos->second;
246  if (table != 0)
247  {
248  sigma = table->FindValue(ekin);
249  }
250  }
251  else
252  {
253  G4Exception("G4MicroElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
254  }
255  }
256 
257  if (verboseLevel > 3)
258  {
259  G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
260  G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
261  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
262  }
263 
264  }
265 
266  return sigma*density;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 
271 void G4MicroElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
272  const G4MaterialCutsCouple* /*couple*/,
273  const G4DynamicParticle* aDynamicElectron,
274  G4double,
275  G4double)
276 {
277 
278  if (verboseLevel > 3)
279  G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
280 
281  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
282 
283  if (electronEnergy0 < killBelowEnergy)
284  {
288  return ;
289  }
290 
291  if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
292  {
293  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
294 
295  G4double phi = 2. * pi * G4UniformRand();
296 
297  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
298  G4ThreeVector xVers = zVers.orthogonal();
299  G4ThreeVector yVers = zVers.cross(xVers);
300 
301  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
302  G4double yDir = xDir;
303  xDir *= std::cos(phi);
304  yDir *= std::sin(phi);
305 
306  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
307 
309 
311  }
312 
313 }
314 
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316 
318  (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
319 {
320 
321  G4double theta = 0.;
322  G4double valueT1 = 0;
323  G4double valueT2 = 0;
324  G4double valueE21 = 0;
325  G4double valueE22 = 0;
326  G4double valueE12 = 0;
327  G4double valueE11 = 0;
328  G4double xs11 = 0;
329  G4double xs12 = 0;
330  G4double xs21 = 0;
331  G4double xs22 = 0;
332 
333 
334  if (particleDefinition == G4Electron::ElectronDefinition())
335  {
336  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
337  std::vector<double>::iterator t1 = t2-1;
338 
339  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
340  std::vector<double>::iterator e11 = e12-1;
341 
342  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
343  std::vector<double>::iterator e21 = e22-1;
344 
345  valueT1 =*t1;
346  valueT2 =*t2;
347  valueE21 =*e21;
348  valueE22 =*e22;
349  valueE12 =*e12;
350  valueE11 =*e11;
351 
352  xs11 = eDiffCrossSectionData[valueT1][valueE11];
353  xs12 = eDiffCrossSectionData[valueT1][valueE12];
354  xs21 = eDiffCrossSectionData[valueT2][valueE21];
355  xs22 = eDiffCrossSectionData[valueT2][valueE22];
356 
357 }
358 
359  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
360 
361  theta = QuadInterpolator( valueE11, valueE12,
362  valueE21, valueE22,
363  xs11, xs12,
364  xs21, xs22,
365  valueT1, valueT2,
366  k, integrDiff );
367 
368  return theta;
369 }
370 
371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
372 
374  G4double e2,
375  G4double e,
376  G4double xs1,
377  G4double xs2)
378 {
379  G4double d1 = std::log(xs1);
380  G4double d2 = std::log(xs2);
381  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
382  return value;
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386 
388  G4double e2,
389  G4double e,
390  G4double xs1,
391  G4double xs2)
392 {
393  G4double d1 = xs1;
394  G4double d2 = xs2;
395  G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
396  return value;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
402  G4double e2,
403  G4double e,
404  G4double xs1,
405  G4double xs2)
406 {
407  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
408  G4double b = std::log10(xs2) - a*std::log10(e2);
409  G4double sigma = a*std::log10(e) + b;
410  G4double value = (std::pow(10.,sigma));
411  return value;
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
415 
417  G4double e21, G4double e22,
418  G4double xs11, G4double xs12,
419  G4double xs21, G4double xs22,
420  G4double t1, G4double t2,
421  G4double t, G4double e)
422 {
423  // Log-Log
424 /*
425  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
426  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
427  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
428 
429 
430  // Lin-Log
431  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
432  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
433  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
434 */
435 
436  // Lin-Lin
437  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
438  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
439  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
440 
441  return value;
442 }
443 
444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
445 
447 {
448  G4double integrdiff=0;
449  G4double uniformRand=G4UniformRand();
450  integrdiff = uniformRand;
451 
452  G4double theta=0.;
453  G4double cosTheta=0.;
454  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
455 
456  cosTheta= std::cos(theta*pi/180);
457 
458  return cosTheta;
459 }
static const double cm
Definition: G4SIunits.hh:106
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:623
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static const double MeV
Definition: G4SIunits.hh:193
static const G4double d1
G4MicroElecElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecElasticModel")
G4double Theta(G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:616
std::vector< double > eTdummyVec
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4double pi
static const G4double e2
G4double a
Definition: TRTMaterials.hh:39
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
G4double density
Definition: TRTMaterials.hh:39
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
const G4ThreeVector & GetMomentumDirection() const
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
static const double eV
Definition: G4SIunits.hh:194
G4ParticleChangeForGamma * fParticleChangeForGamma
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4double RandomizeCosTheta(G4double k)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
static const G4double d2
#define DBL_MAX
Definition: templates.hh:83
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134