Geant4_10
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 
62  nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
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;
68  SetLowEnergyLimit(lowEnergyLimit);
69  SetHighEnergyLimit(highEnergyLimit);
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 
120  if (LowEnergyLimit() < lowEnergyLimit)
121  {
122  G4cout << "G4MuElecElasticModel: low energy limit increased from " <<
123  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
124  SetLowEnergyLimit(lowEnergyLimit);
125  }
126 
127  if (HighEnergyLimit() > highEnergyLimit)
128  {
129  G4cout << "G4MuElecElasticModel: high energy limit decreased from " <<
130  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
131  SetHighEnergyLimit(highEnergyLimit);
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 
318 G4double G4MuElecElasticModel::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 G4MuElecElasticModel::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 = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
383  return value;
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
387 
388 G4double G4MuElecElasticModel::LogLogInterpolate(G4double e1,
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 
403 G4double G4MuElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
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 
419 G4double G4MuElecElasticModel::RandomizeCosTheta(G4double k)
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 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)
virtual G4double FindValue(G4double e, G4int componentId=0) const
tuple a
Definition: test.py:11
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
const char * p
Definition: xmltok.h:285
virtual G4bool LoadData(const G4String &argFileName)
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)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
string material
Definition: eplot.py:19
G4double density
Definition: TRTMaterials.hh:39
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
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)
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
Hep3Vector orthogonal() const
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:231
G4ParticleChangeForGamma * fParticleChangeForGamma
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