Geant4  10.02.p01
G4MuElecElasticModel.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 // G4MuElecElasticModel.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 "G4MuElecElasticModel.hh"
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 
53  G4cout << G4endl;
54  G4cout << "*******************************************************************************" << G4endl;
55  G4cout << "*******************************************************************************" << G4endl;
56  G4cout << " The name of the class G4MuElecElasticModel is changed to G4MicroElecElasticModel. " << G4endl;
57  G4cout << " The obsolete class will be REMOVED with the next release of Geant4. " << G4endl;
58  G4cout << "*******************************************************************************" << G4endl;
59  G4cout << "*******************************************************************************" << G4endl;
60  G4cout << G4endl;
61 
63 
64  killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation
65  lowEnergyLimit = 0 * eV;
66  lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV
67  highEnergyLimit = 100. * MeV;
70 
71  verboseLevel= 0;
72  // Verbosity scale:
73  // 0 = nothing
74  // 1 = warning for energy non-conservation
75  // 2 = details of energy budget
76  // 3 = calculation of cross sections, file openings, sampling of atoms
77  // 4 = entering in methods
78 
79  if( verboseLevel>0 )
80  {
81  G4cout << "MuElec Elastic model is constructed " << G4endl
82  << "Energy range: "
83  << lowEnergyLimit / eV << " eV - "
84  << highEnergyLimit / keV << " keV"
85  << G4endl;
86  }
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
93 {
94  // For total cross section
95 
96  std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
97  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
98  {
99  G4MuElecCrossSectionDataSet* table = pos->second;
100  delete table;
101  }
102 
103  // For final state
104 
105  eVecm.clear();
106 
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
112  const G4DataVector& /*cuts*/)
113 {
114 
115  if (verboseLevel > 3)
116  G4cout << "Calling G4MuElecElasticModel::Initialise()" << G4endl;
117 
118  // Energy limits
119 
121  {
122  G4cout << "G4MuElecElasticModel: low energy limit increased from " <<
123  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
125  }
126 
128  {
129  G4cout << "G4MuElecElasticModel: high energy limit decreased from " <<
130  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
132  }
133 
134  // Reading of data files
135 
136  G4double scaleFactor = 1e-18 * cm * cm;
137 
138  G4String fileElectron("microelec/sigma_elastic_e_Si");
139 
142 
143  // For total cross section
144 
145  electron = electronDef->GetParticleName();
146 
147  tableFile[electron] = fileElectron;
148 
150  tableE->LoadData(fileElectron);
151  tableData[electron] = tableE;
152 
153  // For final state
154 
155  char *path = getenv("G4LEDATA");
156 
157  if (!path)
158  {
159  G4Exception("G4MuElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
160  return;
161  }
162 
163  std::ostringstream eFullFileName;
164  eFullFileName << path << "/microelec/sigmadiff_elastic_e_Si.dat";
165  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
166 
167  if (!eDiffCrossSection)
168  G4Exception("G4MuElecElasticModel::Initialise","em0003",FatalException,"Missing data file: /microelec/sigmadiff_elastic_e_Si.dat");
169 
170  eTdummyVec.push_back(0.);
171 
172  while(!eDiffCrossSection.eof())
173  {
174  double tDummy;
175  double eDummy;
176  eDiffCrossSection>>tDummy>>eDummy;
177 
178  // SI : mandatory eVecm initialization
179  if (tDummy != eTdummyVec.back())
180  {
181  eTdummyVec.push_back(tDummy);
182  eVecm[tDummy].push_back(0.);
183  }
184 
185  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
186 
187  // SI : only if not end of file reached !
188  if (!eDiffCrossSection.eof()) eDiffCrossSectionData[tDummy][eDummy]*=scaleFactor;
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 MuElec Elastic model" << G4endl;
198 
199  if( verboseLevel>0 )
200  {
201  G4cout << "MuElec Elastic model is initialized " << G4endl
202  << "Energy range: "
203  << LowEnergyLimit() / eV << " eV - "
204  << HighEnergyLimit() / keV << " keV"
205  << G4endl;
206  }
207 
208  if (isInitialised) { return; }
210  isInitialised = true;
211 
212  // InitialiseElementSelectors(particle,cuts);
213 }
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216 
218  const G4ParticleDefinition* p,
219  G4double ekin,
220  G4double,
221  G4double)
222 {
223  if (verboseLevel > 3)
224  G4cout << "Calling CrossSectionPerVolume() of G4MuElecElasticModel" << G4endl;
225 
226  // Calculate total cross section for model
227 
228  G4double sigma=0;
229 
231 
232  if (material == nistSi || material->GetBaseMaterial() == nistSi)
233  {
234  const G4String& particleName = p->GetParticleName();
235 
236  if (ekin < highEnergyLimit)
237  {
238  //SI : XS must not be zero otherwise sampling of secondaries method ignored
239  if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
240  //
241 
242  std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
243  pos = tableData.find(particleName);
244 
245  if (pos != tableData.end())
246  {
247  G4MuElecCrossSectionDataSet* table = pos->second;
248  if (table != 0)
249  {
250  sigma = table->FindValue(ekin);
251  }
252  }
253  else
254  {
255  G4Exception("G4MuElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
256  }
257  }
258 
259  if (verboseLevel > 3)
260  {
261  G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
262  G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
263  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
264  }
265 
266  }
267 
268  return sigma*density;
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272 
273 void G4MuElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
274  const G4MaterialCutsCouple* /*couple*/,
275  const G4DynamicParticle* aDynamicElectron,
276  G4double,
277  G4double)
278 {
279 
280  if (verboseLevel > 3)
281  G4cout << "Calling SampleSecondaries() of G4MuElecElasticModel" << G4endl;
282 
283  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
284 
285  if (electronEnergy0 < killBelowEnergy)
286  {
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 = std::exp(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 a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
395  G4double b = std::log10(xs2) - a*std::log10(e2);
396  G4double sigma = a*std::log10(e) + b;
397  G4double value = (std::pow(10.,sigma));
398  return value;
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
402 
404  G4double e21, G4double e22,
405  G4double xs11, G4double xs12,
406  G4double xs21, G4double xs22,
407  G4double t1, G4double t2,
408  G4double t, G4double e)
409 {
410 // Lin-Log
411  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
412  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
413  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
414  return value;
415 }
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
418 
420 {
421  G4double integrdiff=0;
422  G4double uniformRand=G4UniformRand();
423  integrdiff = uniformRand;
424 
425  G4double theta=0.;
426  G4double cosTheta=0.;
427  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
428 
429  cosTheta= std::cos(theta*pi/180);
430 
431  return cosTheta;
432 }
static const double cm
Definition: G4SIunits.hh:118
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
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
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
std::vector< double > eTdummyVec
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
static const G4double e2
virtual G4bool LoadData(const G4String &argFileName)
G4double a
Definition: TRTMaterials.hh:39
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double RandomizeCosTheta(G4double k)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4double density
Definition: TRTMaterials.hh:39
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
const G4ThreeVector & GetMomentumDirection() const
G4MuElecElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MuElecElasticModel")
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double Theta(G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
TriDimensionMap eDiffCrossSectionData
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
static const double eV
Definition: G4SIunits.hh:212
static const double pi
Definition: G4SIunits.hh:74
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4ParticleChangeForGamma * fParticleChangeForGamma
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static const G4double d2
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)
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134