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