Geant4_10
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 {
52  nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
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;
58  SetLowEnergyLimit(lowEnergyLimit);
59  SetHighEnergyLimit(highEnergyLimit);
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 / keV << " keV"
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 
110  if (LowEnergyLimit() < lowEnergyLimit)
111  {
112  G4cout << "G4MicroElecElasticModel: low energy limit increased from " <<
113  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
114  SetLowEnergyLimit(lowEnergyLimit);
115  }
116 
117  if (HighEnergyLimit() > highEnergyLimit)
118  {
119  G4cout << "G4MicroElecElasticModel: high energy limit decreased from " <<
120  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
121  SetHighEnergyLimit(highEnergyLimit);
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_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_elastic_e_Si.dat");
159 
160  eTdummyVec.push_back(0.);
161 
162  while(!eDiffCrossSection.eof())
163  {
164  double tDummy;
165  double eDummy;
166  eDiffCrossSection>>tDummy>>eDummy;
167 
168  // SI : mandatory eVecm initialization
169  if (tDummy != eTdummyVec.back())
170  {
171  eTdummyVec.push_back(tDummy);
172  eVecm[tDummy].push_back(0.);
173  }
174 
175  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
176 
177  // SI : only if not end of file reached !
178  if (!eDiffCrossSection.eof()) eDiffCrossSectionData[tDummy][eDummy]*=scaleFactor;
179 
180  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
181 
182  }
183 
184  // End final state
185 
186  if (verboseLevel > 2)
187  G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
188 
189  if( verboseLevel>0 )
190  {
191  G4cout << "MicroElec Elastic model is initialized " << G4endl
192  << "Energy range: "
193  << LowEnergyLimit() / eV << " eV - "
194  << HighEnergyLimit() / keV << " keV"
195  << G4endl;
196  }
197 
198  if (isInitialised) { return; }
200  isInitialised = true;
201 
202  // InitialiseElementSelectors(particle,cuts);
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 
208  const G4ParticleDefinition* p,
209  G4double ekin,
210  G4double,
211  G4double)
212 {
213  if (verboseLevel > 3)
214  G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
215 
216  // Calculate total cross section for model
217 
218  G4double sigma=0;
219 
221 
222  if (material == nistSi || material->GetBaseMaterial() == nistSi)
223  {
224  const G4String& particleName = p->GetParticleName();
225 
226  if (ekin < highEnergyLimit)
227  {
228  //SI : XS must not be zero otherwise sampling of secondaries method ignored
229  if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
230  //
231 
232  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
233  pos = tableData.find(particleName);
234 
235  if (pos != tableData.end())
236  {
237  G4MicroElecCrossSectionDataSet* table = pos->second;
238  if (table != 0)
239  {
240  sigma = table->FindValue(ekin);
241  }
242  }
243  else
244  {
245  G4Exception("G4MicroElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
246  }
247  }
248 
249  if (verboseLevel > 3)
250  {
251  G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
252  G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
253  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
254  }
255 
256  }
257 
258  return sigma*density;
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 
263 void G4MicroElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
264  const G4MaterialCutsCouple* /*couple*/,
265  const G4DynamicParticle* aDynamicElectron,
266  G4double,
267  G4double)
268 {
269 
270  if (verboseLevel > 3)
271  G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
272 
273  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
274 
275  if (electronEnergy0 < killBelowEnergy)
276  {
279  return ;
280  }
281 
282  if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
283  {
284  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
285 
286  G4double phi = 2. * pi * G4UniformRand();
287 
288  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
289  G4ThreeVector xVers = zVers.orthogonal();
290  G4ThreeVector yVers = zVers.cross(xVers);
291 
292  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
293  G4double yDir = xDir;
294  xDir *= std::cos(phi);
295  yDir *= std::sin(phi);
296 
297  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
298 
300 
302  }
303 
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
308 G4double G4MicroElecElasticModel::Theta
309  (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
310 {
311 
312  G4double theta = 0.;
313  G4double valueT1 = 0;
314  G4double valueT2 = 0;
315  G4double valueE21 = 0;
316  G4double valueE22 = 0;
317  G4double valueE12 = 0;
318  G4double valueE11 = 0;
319  G4double xs11 = 0;
320  G4double xs12 = 0;
321  G4double xs21 = 0;
322  G4double xs22 = 0;
323 
324 
325  if (particleDefinition == G4Electron::ElectronDefinition())
326  {
327  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
328  std::vector<double>::iterator t1 = t2-1;
329 
330  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
331  std::vector<double>::iterator e11 = e12-1;
332 
333  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
334  std::vector<double>::iterator e21 = e22-1;
335 
336  valueT1 =*t1;
337  valueT2 =*t2;
338  valueE21 =*e21;
339  valueE22 =*e22;
340  valueE12 =*e12;
341  valueE11 =*e11;
342 
343  xs11 = eDiffCrossSectionData[valueT1][valueE11];
344  xs12 = eDiffCrossSectionData[valueT1][valueE12];
345  xs21 = eDiffCrossSectionData[valueT2][valueE21];
346  xs22 = eDiffCrossSectionData[valueT2][valueE22];
347 
348 }
349 
350  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
351 
352  theta = QuadInterpolator( valueE11, valueE12,
353  valueE21, valueE22,
354  xs11, xs12,
355  xs21, xs22,
356  valueT1, valueT2,
357  k, integrDiff );
358 
359  return theta;
360 }
361 
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
363 
364 G4double G4MicroElecElasticModel::LinLogInterpolate(G4double e1,
365  G4double e2,
366  G4double e,
367  G4double xs1,
368  G4double xs2)
369 {
370  G4double d1 = std::log(xs1);
371  G4double d2 = std::log(xs2);
372  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
373  return value;
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
377 
378 G4double G4MicroElecElasticModel::LogLogInterpolate(G4double e1,
379  G4double e2,
380  G4double e,
381  G4double xs1,
382  G4double xs2)
383 {
384  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
385  G4double b = std::log10(xs2) - a*std::log10(e2);
386  G4double sigma = a*std::log10(e) + b;
387  G4double value = (std::pow(10.,sigma));
388  return value;
389 }
390 
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392 
393 G4double G4MicroElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
394  G4double e21, G4double e22,
395  G4double xs11, G4double xs12,
396  G4double xs21, G4double xs22,
397  G4double t1, G4double t2,
398  G4double t, G4double e)
399 {
400 // Lin-Log
401  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
402  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
403  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
404  return value;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
409 G4double G4MicroElecElasticModel::RandomizeCosTheta(G4double k)
410 {
411  G4double integrdiff=0;
412  G4double uniformRand=G4UniformRand();
413  integrdiff = uniformRand;
414 
415  G4double theta=0.;
416  G4double cosTheta=0.;
417  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
418 
419  cosTheta= std::cos(theta*pi/180);
420 
421  return cosTheta;
422 }
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
TTree * t1
Definition: plottest35.C:26
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4MicroElecElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecElasticModel")
tuple a
Definition: test.py:11
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
const char * p
Definition: xmltok.h:285
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
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:683
string material
Definition: eplot.py:19
G4double density
Definition: TRTMaterials.hh:39
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
const G4ThreeVector & GetMomentumDirection() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
Hep3Vector unit() const
G4ParticleChangeForGamma * fParticleChangeForGamma
Hep3Vector orthogonal() const
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:231
const XML_Char int const XML_Char * value
Definition: expat.h:331
TTree * t2
Definition: plottest35.C:36
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121