ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MicroElecElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MicroElecElasticModel.cc
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 
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 
102 void G4MicroElecElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
103  const G4DataVector& /*cuts*/)
104 {
105 
106  if (verboseLevel > 3)
107  G4cout << "Calling G4MicroElecElasticModel::Initialise()" << G4endl;
108 
109  // Energy limits
110 
111  if (LowEnergyLimit() < lowEnergyLimit)
112  {
113  G4cout << "G4MicroElecElasticModel: low energy limit increased from " <<
114  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
115  SetLowEnergyLimit(lowEnergyLimit);
116  }
117 
118  if (HighEnergyLimit() > highEnergyLimit)
119  {
120  G4cout << "G4MicroElecElasticModel: high energy limit decreased from " <<
121  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
122  SetHighEnergyLimit(highEnergyLimit);
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 
131  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
132  G4String electron;
133 
134  // For total cross section
135 
136  electron = electronDef->GetParticleName();
137 
138  tableFile[electron] = fileElectron;
139 
140  G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141  tableE->LoadData(fileElectron);
142  tableData[electron] = tableE;
143 
144  // For final state
145 
146  char *path = std::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; }
209  fParticleChangeForGamma = GetParticleChangeForGamma();
210  isInitialised = true;
211 
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215 
216 G4double G4MicroElecElasticModel::CrossSectionPerVolume(const G4Material* material,
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 
229  G4double density = material->GetTotNbOfAtomsPerVolume();
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  {
286  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
287  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
288  fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
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 
309  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
310 
311  fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
312  }
313 
314 }
315 
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317 
318 G4double G4MicroElecElasticModel::Theta
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 
374 G4double G4MicroElecElasticModel::LinLogInterpolate(G4double e1,
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 
388 G4double G4MicroElecElasticModel::LinLinInterpolate(G4double e1,
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 
402 G4double G4MicroElecElasticModel::LogLogInterpolate(G4double e1,
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 
417 G4double G4MicroElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
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 
447 G4double G4MicroElecElasticModel::RandomizeCosTheta(G4double k)
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 }