Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 #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 {
53 
54  G4cout << G4endl;
55  G4cout << "*******************************************************************************" << G4endl;
56  G4cout << "*******************************************************************************" << G4endl;
57  G4cout << " The name of the class G4MuElecElasticModel is changed to G4MicroElecElasticModel. " << G4endl;
58  G4cout << " The obsolete class will be REMOVED with the next release of Geant4. " << G4endl;
59  G4cout << "*******************************************************************************" << G4endl;
60  G4cout << "*******************************************************************************" << G4endl;
61  G4cout << G4endl;
62 
63  nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
64 
65  killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation
66  lowEnergyLimit = 0 * eV;
67  lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV
68  highEnergyLimit = 100. * MeV;
69  SetLowEnergyLimit(lowEnergyLimit);
70  SetHighEnergyLimit(highEnergyLimit);
71 
72  verboseLevel= 0;
73  // Verbosity scale:
74  // 0 = nothing
75  // 1 = warning for energy non-conservation
76  // 2 = details of energy budget
77  // 3 = calculation of cross sections, file openings, sampling of atoms
78  // 4 = entering in methods
79 
80  if( verboseLevel>0 )
81  {
82  G4cout << "MuElec Elastic model is constructed " << G4endl
83  << "Energy range: "
84  << lowEnergyLimit / eV << " eV - "
85  << highEnergyLimit / keV << " keV"
86  << G4endl;
87  }
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94 {
95  // For total cross section
96 
97  std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
98  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99  {
100  G4MuElecCrossSectionDataSet* table = pos->second;
101  delete table;
102  }
103 
104  // For final state
105 
106  eVecm.clear();
107 
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
113  const G4DataVector& /*cuts*/)
114 {
115 
116  if (verboseLevel > 3)
117  G4cout << "Calling G4MuElecElasticModel::Initialise()" << G4endl;
118 
119  // Energy limits
120 
121  if (LowEnergyLimit() < lowEnergyLimit)
122  {
123  G4cout << "G4MuElecElasticModel: low energy limit increased from " <<
124  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
125  SetLowEnergyLimit(lowEnergyLimit);
126  }
127 
128  if (HighEnergyLimit() > highEnergyLimit)
129  {
130  G4cout << "G4MuElecElasticModel: high energy limit decreased from " <<
131  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
132  SetHighEnergyLimit(highEnergyLimit);
133  }
134 
135  // Reading of data files
136 
137  G4double scaleFactor = 1e-18 * cm * cm;
138 
139  G4String fileElectron("microelec/sigma_elastic_e_Si");
140 
143 
144  // For total cross section
145 
146  electron = electronDef->GetParticleName();
147 
148  tableFile[electron] = fileElectron;
149 
151  tableE->LoadData(fileElectron);
152  tableData[electron] = tableE;
153 
154  // For final state
155 
156  char *path = getenv("G4LEDATA");
157 
158  if (!path)
159  {
160  G4Exception("G4MuElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
161  return;
162  }
163 
164  std::ostringstream eFullFileName;
165  eFullFileName << path << "/microelec/sigmadiff_elastic_e_Si.dat";
166  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
167 
168  if (!eDiffCrossSection)
169  G4Exception("G4MuElecElasticModel::Initialise","em0003",FatalException,"Missing data file: /microelec/sigmadiff_elastic_e_Si.dat");
170 
171  eTdummyVec.push_back(0.);
172 
173  while(!eDiffCrossSection.eof())
174  {
175  double tDummy;
176  double eDummy;
177  eDiffCrossSection>>tDummy>>eDummy;
178 
179  // SI : mandatory eVecm initialization
180  if (tDummy != eTdummyVec.back())
181  {
182  eTdummyVec.push_back(tDummy);
183  eVecm[tDummy].push_back(0.);
184  }
185 
186  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
187 
188  // SI : only if not end of file reached !
189  if (!eDiffCrossSection.eof()) eDiffCrossSectionData[tDummy][eDummy]*=scaleFactor;
190 
191  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
192 
193  }
194 
195  // End final state
196 
197  if (verboseLevel > 2)
198  G4cout << "Loaded cross section files for MuElec Elastic model" << G4endl;
199 
200  if( verboseLevel>0 )
201  {
202  G4cout << "MuElec Elastic model is initialized " << G4endl
203  << "Energy range: "
204  << LowEnergyLimit() / eV << " eV - "
205  << HighEnergyLimit() / keV << " keV"
206  << G4endl;
207  }
208 
209  if (isInitialised) { return; }
211  isInitialised = true;
212 
213  // InitialiseElementSelectors(particle,cuts);
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217 
219  const G4ParticleDefinition* p,
220  G4double ekin,
221  G4double,
222  G4double)
223 {
224  if (verboseLevel > 3)
225  G4cout << "Calling CrossSectionPerVolume() of G4MuElecElasticModel" << G4endl;
226 
227  // Calculate total cross section for model
228 
229  G4double sigma=0;
230 
231  G4double density = material->GetTotNbOfAtomsPerVolume();
232 
233  if (material == nistSi || material->GetBaseMaterial() == nistSi)
234  {
235  const G4String& particleName = p->GetParticleName();
236 
237  if (ekin < highEnergyLimit)
238  {
239  //SI : XS must not be zero otherwise sampling of secondaries method ignored
240  if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
241  //
242 
243  std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
244  pos = tableData.find(particleName);
245 
246  if (pos != tableData.end())
247  {
248  G4MuElecCrossSectionDataSet* table = pos->second;
249  if (table != 0)
250  {
251  sigma = table->FindValue(ekin);
252  }
253  }
254  else
255  {
256  G4Exception("G4MuElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
257  }
258  }
259 
260  if (verboseLevel > 3)
261  {
262  G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
263  G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl;
264  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl;
265  }
266 
267  }
268 
269  return sigma*density;
270 }
271 
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273 
274 void G4MuElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
275  const G4MaterialCutsCouple* /*couple*/,
276  const G4DynamicParticle* aDynamicElectron,
277  G4double,
278  G4double)
279 {
280 
281  if (verboseLevel > 3)
282  G4cout << "Calling SampleSecondaries() of G4MuElecElasticModel" << G4endl;
283 
284  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
285 
286  if (electronEnergy0 < killBelowEnergy)
287  {
290  return ;
291  }
292 
293  if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
294  {
295  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
296 
297  G4double phi = 2. * pi * G4UniformRand();
298 
299  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
300  G4ThreeVector xVers = zVers.orthogonal();
301  G4ThreeVector yVers = zVers.cross(xVers);
302 
303  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
304  G4double yDir = xDir;
305  xDir *= std::cos(phi);
306  yDir *= std::sin(phi);
307 
308  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
309 
311 
313  }
314 
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
318 
319 G4double G4MuElecElasticModel::Theta
320  (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
321 {
322 
323  G4double theta = 0.;
324  G4double valueT1 = 0;
325  G4double valueT2 = 0;
326  G4double valueE21 = 0;
327  G4double valueE22 = 0;
328  G4double valueE12 = 0;
329  G4double valueE11 = 0;
330  G4double xs11 = 0;
331  G4double xs12 = 0;
332  G4double xs21 = 0;
333  G4double xs22 = 0;
334 
335 
336  if (particleDefinition == G4Electron::ElectronDefinition())
337  {
338  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
339  std::vector<double>::iterator t1 = t2-1;
340 
341  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
342  std::vector<double>::iterator e11 = e12-1;
343 
344  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
345  std::vector<double>::iterator e21 = e22-1;
346 
347  valueT1 =*t1;
348  valueT2 =*t2;
349  valueE21 =*e21;
350  valueE22 =*e22;
351  valueE12 =*e12;
352  valueE11 =*e11;
353 
354  xs11 = eDiffCrossSectionData[valueT1][valueE11];
355  xs12 = eDiffCrossSectionData[valueT1][valueE12];
356  xs21 = eDiffCrossSectionData[valueT2][valueE21];
357  xs22 = eDiffCrossSectionData[valueT2][valueE22];
358 
359 }
360 
361  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
362 
363  theta = QuadInterpolator( valueE11, valueE12,
364  valueE21, valueE22,
365  xs11, xs12,
366  xs21, xs22,
367  valueT1, valueT2,
368  k, integrDiff );
369 
370  return theta;
371 }
372 
373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374 
375 G4double G4MuElecElasticModel::LinLogInterpolate(G4double e1,
376  G4double e2,
377  G4double e,
378  G4double xs1,
379  G4double xs2)
380 {
381  G4double d1 = std::log(xs1);
382  G4double d2 = std::log(xs2);
383  G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
384  return value;
385 }
386 
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388 
389 G4double G4MuElecElasticModel::LogLogInterpolate(G4double e1,
390  G4double e2,
391  G4double e,
392  G4double xs1,
393  G4double xs2)
394 {
395  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
396  G4double b = std::log10(xs2) - a*std::log10(e2);
397  G4double sigma = a*std::log10(e) + b;
398  G4double value = (std::pow(10.,sigma));
399  return value;
400 }
401 
402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
403 
404 G4double G4MuElecElasticModel::QuadInterpolator(G4double e11, G4double e12,
405  G4double e21, G4double e22,
406  G4double xs11, G4double xs12,
407  G4double xs21, G4double xs22,
408  G4double t1, G4double t2,
409  G4double t, G4double e)
410 {
411 // Lin-Log
412  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
413  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
414  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
415  return value;
416 }
417 
418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
419 
420 G4double G4MuElecElasticModel::RandomizeCosTheta(G4double k)
421 {
422  G4double integrdiff=0;
423  G4double uniformRand=G4UniformRand();
424  integrdiff = uniformRand;
425 
426  G4double theta=0.;
427  G4double cosTheta=0.;
428  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
429 
430  cosTheta= std::cos(theta*pi/180);
431 
432  return cosTheta;
433 }
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
static const G4double d2
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:732
string material
Definition: eplot.py:19
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:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double eV
Definition: G4SIunits.hh:215
G4MuElecElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MuElecElasticModel")
static const G4double d1
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 G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
Hep3Vector unit() const
tuple t1
Definition: plottest35.py:33
Hep3Vector orthogonal() const
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4ParticleChangeForGamma * fParticleChangeForGamma
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
static constexpr double keV
Definition: G4SIunits.hh:216
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132